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This  research  program  investigated  a  new  tool  for  the  analysis  of  elec¬ 
tromagnetic  coupling  and  shielding  problems:  the  finite-difference, 
time-domain  (FD-TD)  solution  of  Maxwell's  equations.  The  objective  of 
the  program  was  to  evaluate  the  suitability  of  the  FD-TD  method  to 
determine  the  amount  of  electromagnetic  coupling  through  an  aperture  into 
an  enclosed  conducting  container  and  the  interaction  and  coupling  of  the 
penetrating  fields  with  internal  electronics.  Two  specific  container 
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models  were  used  for  the  evaluation.  The  first,  a  conducting  cylinder 
with  one  open  end;  the  other,  the  guidance  section  of  a  missile.  Each 
of  these  two  configurations  was  modeled  to  calculate  the  electromagnetic 
field  coupled  into  the  structure. 

The  following  specific  questions  were  addressed  during  this  program: 

1.  Can  the  FD-TD  method  accurately  model  electromagnetic  coupling 
into  a  conducting  structure  for  arbitrary  angles  of  incidence  and 
arbitrary  wave  polarization? 

2.  Can  the  FD-TD  method  accurately  model  a  complex  structure  with 
both  irregular-shaped  apertures  and  dielectric  or  permeable  materials 
within  the  interior? 

3.  Can  the  FD-TD  method  predict  the  voltage,  current,  or  power  that 
the  coupled  field  may  indice  on  wires  or  cable  bundles  inside  of  a  metal 
structure? 

4.  What  are  the  size  limits  of  an  object  that  can  be  modeled  using 
the  FD-TD  method? 

5.  What  is  the  capability  of  the  FD-TD  method  to  interface  with 
other  analysis  techniques,  such  as  "Method  of  Moments? 

Overall,  this  program  showed  that  the  FD-TD  method  can  be  success¬ 
fully  applied  to  electromagnetic  coupling  problems  involving  conducting 
structures  with  hole  and  sleeve-type  apertures  illuminated  by  a  plane 
wave  having  an  arbitrary  polarization  and  angle  of  incidence.  Further, 
the  FD-TD  method  can  be  applied  to  complex  cavity-like  structures  having 
internal  metal  and  dielectric  materials,  as  well  as  connecting  wires. 
Accuracy  of  the  FD-TD  results  was  very  good  relative  to  the  uncertainties 
of  available  experimental  and  numerical-theory  approaches.  Convergence 
of  the  electromagnetic  fields  to  the  sinusoidal  steady  state  occurred 
within  about  3  cycles  of  the  incident  wave  when  a  slight  value  of  iso¬ 
tropic  loss  was  assigned  to  the  interior  of  the  structures  modeled.  This 
resulted  in  program  central  processor  times  of  less  than  5  minutes  for 
FD-TD  lattices  containing  as  many  as  1.5  million  unknowns  time-stepped  to 
nmax  =  using  the  Control  Data  STAR-100  and  Cyber  203  computers. 

This  program  also  established  the  feasibility  of  a  hybrid  MOM/FD-TD 
analysis  technique  based  upon  the  use  of  Schelkunoff 's  equivalent 
electric  current  theorem  at  apertures  of  a  cavity.  This  hybrid  technique 
was  shown  to  give  consistent  results  for  structures  illuminated  at  both 
oblique  and  axial  incidence  by  plane  waves.  Further,  this  technique  was 
shown  to  give  consistent  results  for  apertures  strongly  coupled  to 
internal  wires  or  other  metal  and  dielectric  structures. 

This  program  has  shown  that  the  pure  FD-TD  and  hybrid  MOM/FD-TD 
methods  have  great  promise  for  applications  involving  complex  conducting 
and  dielectric  structures  illuminated  by  plane  waves  at  arbitrary  angles 
of  incidence  and  polarization.  Further,  great  promise  is  shown  for  those 
structures  that  are  simultaneously  electrically  large  compared  to  a  wave¬ 
length  and  penetrated  by  locally  complex  cavity-backed  apertures  having 
dimensions  comparable  to  a  wavelength.  Such  structures  requiring 
resolution  of  both  large  and  small  details  at  the  same  time  have  not  been 
well  treated  by  any  one  previous  analytical  or  numerical  approach. 
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1.0  BACKGROUND 


Electromagnetic  penetration  problems  are  difficult  to  treat  with  many 
analytical  or  numerical  methods  because  of  the  inability  of  these  methods  to 
simply  deal  with  the  effects  of  structure  apertures,  curvatures,  corners, 
and  internal  contents.  Usually,  only  relatively  simple  geometries  and  aper¬ 
tures  are  studied  in  an  attempt  to  gain  insight  into  the  key  penetration 
mechanisms  and  to  allow  an  indirect  estimate  of  the  penetration  for  more 
complicated  problems. 


In  RADC  Contracts  F30602-77-C-0163  and  F30602-79-C-0039,  I  IT  Research 
Institute  (IITRI)  investigated  the  application  of  a  new  approach  for  the 
direct  modeling  of  electromagnetic  interaction  problems:  the  finite-differ¬ 
ence,  time-domain  (FD-TD)  solution  of  Maxwell's  equations.  The  FD-TD  method 
treats  the  irradiation  of  a  structure  as  an  initial  value  problem.  At  t  =  0, 
a  plane  wave  source  of  frequency,  f,  is  assumed  to  be  turned  on.  The  propa¬ 
gation  of  waves  from  this  source  is  simulated  by  solving  a  finite-difference 
analog  of  the  time-dependent  Maxwell's  equations  on  a  lattice  of  cells, 
including  the  structure.  Time-stepping  is  continued  until  the  sinusoidal 
steady  state  is  achieved  at  each  cell.  The  field  envelope,  or  maximum  abso¬ 
lute  value,  during  the  final  half  wave-cycle  of  time-stepping  is  taken  as 
the  magnitude  of  the  phasor  of  the  steady-state  field. 


This  method  has  two  key  advantages  relative  to  available  modeling 
approaches.  First,  it  is  simple  to  implement  for  complicated  metal/dielectric 
structures  because  arbitrary  electrical  parameters  can  be  assigned  to  each 
lattice  cell  using  a  data  card  deck.  Second,  its  computer  memory  and  running 
time  requirement  is  not  prohibitive  for  many  complex  structures  of  interest. 

In  the  RADC  work,  IITRI  has  shown  the  FD-TD  method  to  be  capable  of  accurately 
solving  for  hundred  of  thousands  of  unknown  field  components  within  a  few 
minutes  on  an  array-processing  computer.  Consistently,  a  ±  1  dB  accuracy 
relative  to  known  analytical  and  experimental  bench  marks  has  been  achieved 


for  a  variety  of  dielectric  and  metal  geometries. 
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The  objective  of  IITRI's  RAD C  study  was  to  evaluate  the  suitability  of 
the  FD-TD  method  to  determine  the  amount  of  electromagnetic  coupling  through 
an  aperture  into  an  enclosed  conducting  container  and  the  interaction  and 
coupling  of  the  penetrating  fields  with  internal  electronics.  Two  specific 
container  models  were  used  for  the  evaluation.  The  first,  a  conducting  cyl¬ 
inder  with  one  open  end.  The  other,  the  guidance  section  of  a  missile.  Each 
of  these  two  configurations  was  modeled  to  calculate  the  electromagnetic  field 
coupled  into  the  structure. 

The  ultimate  aim  of  research  in  this  area  is  two-fold.  First,  develop 
an  easily  used  general  code  solving  for  the  fields  within  arbitrary  metal/ 
dielectric  structures  having  multiple  apertures,  while  assuring  a  known  bound¬ 
ed  level  of  uncertainty.  And  second,  develop  a  more  sophisticated  intuitive 
understanding  of  basic  wave  penetration  mechanisms  in  time  domain,  such  as 
transient  propagation  through  beyond-cutoff  cavity  interiors,  field  build-up 
at  edges,  and  convergence  to  the  sinusoidal  steady  state. 

2.0  OBJECTIVES  OF  RADC  CONTRACT  F30602-79-C-0039 

The  first-stage  evaluation  of  the  FD-TD  method,  RADC  Contract  F30602-77- 
C-0163,  demonstrated  that  the  FD-TD  method  can  be  successfully  applied  to 
axial  -  incidence,  electromagnetic  wave  coupling  problems  involving  highly 
conducting  structures  with  hole  and  sleeve-type  apertures.  Accuracy  of  the 
FD-TD  results  was  very  good  relative  to  the  uncertainties  of  available  exper¬ 
imental  and  numerical -theory  approaches.  Convergence  of  the  EM  fields  to  the 
sinusoidal  steady  state  occurred  within  about  2  cycles  of  the  incident  wave. 
This  resulted  in  program  running  times  of  3.5  minutes  or  less  on  the  CDC 
STAR-100  for  105  -  cell  lattices. 

The  following  questions  concerning  the  FD-TD  method  remained  unanswered 
at  the  conclusion  of  RADC  Contract  F30602-77-C-01 63 : 

1.  Can  the  FD-TD  method  accurately  model  EM  coupling  into  a  conducting 
structure  for  arbitrary  angles  of  incidence  and  arbitrary  wave 
polarization? 
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2.  Can  the  FD-TD  method  accurately  model  a  complex  structure  with 
both  irregular-shaped  apertures  and  dielectric  or  permeable 
materials  within  the  interior? 

3.  Can  the  FD-TD  method  predict  the  voltage,  current,  or  power 
that  the  coupled  field  may  induce  on  wires  or  cable  bundles 
inside  of  a  metal  structure? 

4.  What  are  the  size  limits  of  an  object  that  can  be  modeled  using 
the  FD-TD  method? 

5.  What  is  the  capability  of  the  FD-TD  method  to  interface  with 
other  analysis  techniques,  such  as  "Method  of  Moments"? 

The  objective  of  the  second-stage  evaluation  of  the  FD-TD  method,  RADC 
Contract  F30602-79-C-0039,  was  to  systematically  address  the  five  questions 
raised  above.  It  was  felt  that  favorable  conclusions  on  most  or  all  of  these 
questions  would  probably  lead  to  eventual  application  of  the  FD-TD  method  to 
complex  EM  coupling  problems  of  interest  to  RADC  and  the  electromagnetic  com¬ 
patibility  community. 

3.0  DESCRIPTION  OF  COUPLING  ANALYSIS  PROBLEMS  CONSIDERED 

3.1  Task  1:  Prediction  of  the  Coupling  into  an  Open-Ended 
_ Cylinder  for  Broadside  Incidence _ 

IITRI  employed  the  FD-TD  technique  to  solve  the  following  electromagnetic 
coupling  problems: 

Interacting  structure  -  Circular  (19.0  cm  diameter),  68.5  cm  long, 
open-ended  aluminum  cylinder  (shown  in  Figures  3  and  4  of  the  Final 
Report) ; 

Incident  wave  -  300  MHz  plane  wave  propagating  perpendicular  to  the 
cylinder  axis  (i.e.,  broadside  to  the  cylinder).  Wave  polarization 
cases  investigated: 

Case  1  -  Transverse  electric  (TE),  i.e.,  wave  electric  field 
perpendicular  to  the  plane  formed  by  the  cylinder 
axis  and  the  incident  wavevector; 
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Case  2  -  Transverse  magnetic  (TM),  i.e.,  wave  magnetic  field 
perpendicular  to  the  plane  formed  by  the  cylinder 
axis  and  the  incident  wavevector; 

Desired  fields  -  Each  component  of  total  E  and  total  H  in  the  axial 
cross-section  plane  of  the  cylinder  down  to  40  cm  from  the  open  end. 
First,  with  the  cross-section  plane  parallel  to  the  incident  trans¬ 
verse  field  component,  and  again  with  the  plane  parallel  to  the 
incident  wavevector; 

Resolution  -  1/3  cm  uniformly  throughout  the  mapping  planes; 

Plotted  values  -  In  decibels  relative  to  an  incident  F  of  1  volt/meter 
and  an  incident  ff  of  1/377  ampere/meter. 

To  solve  these  coupling  problems,  IITRI  suitably  modified  an  existing 
FD-TD  computer  code  to  model  the  cylinder  at  broadside  incidence  for  either 
the  TE  or  TM  cases  of  wave  polarization.  Details  of  the  models  are  given  in 
Sections  3.1  and  3.2  of  the  Final  Report. 

3.2  Task  2;  Prediction  of  the  Coupling  into  a  Loaded 

_ Missile  Guidance  Section _ 

IITRI  employed  the  FD-TD  technique  to  solve  the  following  electromagnetic 
coupling  problems: 

Interacting  structure  -  Missile  guidance  section  with  two  apertures: 
a  circular  one  in  the  nose,  and  a  sleeve  fitting  located  23-1/3  cm  aft. 
Missile  body  geometry  beyond  sleeve  fitting  assumed  to  continue  to  in¬ 
finity  with  constant  cross-section  shape.  Models  of  interior  contents 
investigated: 

Case  1  -  Dielectric  materials  of  the  nose  cone  shell,  including 
separate  specifications  of  the  fiberglas  shell  material, 
phenolic  cup,  and  magnesium  fluoride  infrared  dome  (as 
shown  in  Figures  19  -  21  of  the  Final  Report); 

Case  2  -  All  dielectric  material  of  Case  1  with  the  addition  of 
the  following  metal  and  dielectric  interior  components 
(located  as  shown  in  Figures  24  -  26  of  the  Final  Report): 
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1)  Head  coil  assembly  (assumed  so’,  id  metal); 

2)  Cooled  detector  unit  -  CDU  (assumed  solid  metal); 

3)  Phenolic  ring  around  the  CDU; 

4)  Preamp  can  (metal ) ; 

5)  Wire  connecting  the  CDU  to  the  preamp  can; 

6)  Wire  connecting  the  preamp  can  to  the  metal  backplane; 
and 

7)  Longitudinal  metal  support  rods. 

Incident  wave  -  300  MHz  plane  wave  propagating  down  the  axis  of  the 
structure  toward  its  nose  aperture; 

Desired  fields  -  Each  component  of  total  F  and  total  FT  in  two  axial 
cross-section  planes  of  the  structure.  First,  with  the  cross-section 
plane  parallel  to  the  incident  E\  and  again  with  the  plane  parallel 
to  the  incident  H; 

Resolution  -  1/3  cm  uniformly  throughout  the  mapping  planes; 

Plotted  values  -  In  decibels  relative  to  an  incident  E  of  1  volt/meter 

and  an  incident  H  of  1/377  ampere/meter. 

To  solve  these  coupling  problems,  IITRI  suitably  modified  an  existing 
FD-TD  computer  code  and  composed  two  new  data  card  decks  describing  the  sep¬ 
arate  cases  of  the  internal  structure  of  the  guidance  section.  Details  of 

the  models  are  given  in  Sections  4.1  and  4.2  of  the  Final  Report. 

3.3  Task  3:  Prediction  of  the  Coupling  to  Wires  or  Cable  Bundles 

IITRI  conducted  a  study  of  the  ability  of  the  FD-TD  method  to  predict 
electromagnetic  coupl  ing  to  a  wire  or  cable  bundle  using  Case  2  of  the  missile 
guidance  section  model.  As  stated,  this  case  included  two  connecting  wires; 
1)  Wire  1  -  between  the  cooled  detector  unit  and  the  preamp  can;  and  2) 

Wire  2  -  between  the  preamp  can  and  the  metal  backplane.  Wire  1  and  Wire  2 

were  really  idealizations  of  two  multiconductor  wire  bundles  extended  between 
the  structures  mentioned.  Some  wires  of  each  bundle  terminated  in  electronic 
components;  some  were  small -diameter  cables  with  grounded  shields.  The  goal 
was  to  model  the  direct  metal-to-metal  current  path  through  the  missile  guid¬ 
ance  section.  Therefore,  using  a  simple  single-wire  model  for  each  bundle, 
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only  the  ground-path  current  and  equivalent  common-mode  current  of  each 
bundle  was  modeled  using  the  FD-TD  method. 

Modeling  of  Wire  1  and  Wire  2  was  done  using  the  FD-TD  method  in  a  direct 
manner  by  specifying  each  wire  model  as  a  filamentary  conducting  path  of  real¬ 
istic  location  relative  to  the  other  components  of  the  model.  In  this  manner, 
certain  transmission  line  properties  of  each  wire,  such  as  capacitance  and 
external  inductance,  that  arise  due  to  its  location  relative  to  other  struc¬ 
tures  are  automatical ly  accounted  for  in  the  time-domain  solution  of  Maxwell's 
equations  for  the  whole  model,  which  includes  the  interaction  of  the  wire  and 
all  adjacent  structures.  The  transmission-line  properties  of  each  wire  that 
arise  due  to  internal  impedance  are  accounted  for  by  adjusting  the  model  wire 
conductivity  so  that  these  properties  match  those  of  the  actual  conductor. 

This  procedure  is  discussed  in  Section  5  of  the  Final  Report. 

3.4  Task  4:  Determination  of  the  Size  Limits  of  an  Object 

_ that  can  be  Modeled  by  the  FD-TD  Method _ 

IITRI  conducted  a  study  of  the  maximum  size  of  an  object  that  can  be 
modeled  using  the  FD-TD  method,  given  limitations  of  the  available  computer 
system.  In  general ,  the  computer  resources  needed  by  the  FD-TD  method  in¬ 
crease  for  larger  structures  and  finer  spatial  resolutions,  since  more  memory 
and  processing  time  is  required  to  account  for  higher  levels  of  geometric 
detail.  Factors  such  as  computer  central  memory  size,  data  transmission  rate 
of  on-line  mass  storage  memory,  floating-point  operation  speed,  and  computer 
architecture  were  considered  in  this  study.  Details  of  the  study  are  report¬ 
ed  in  Section  6  of  the  Final  Report. 

3.5  Task  5:  Investigation  of  Interfacing  the  FD-TD  Method 

_ with  Other  Analysis  Techniques _ 

Under  this  task,  IITRI  developed  a  new  hybrid  method-of-moments  (MOM)  / 
FD-TD  technique  based  on  a  novel  use  of  an  aperture  electric  current  field 
equivalence  theorem  due  to  Schelkunoff.  This  hybrid  technique  and  its  appli¬ 
cations,  described  in  Section  7  of  the  Final  Report,  is  basically  an  equiv¬ 
alent  aperture  excitation  method.  This  allows  analysis  of  the  coupling  prob¬ 
lem  in  two  distinct  steps: 
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Step  1  -  Analysis  of  the  relatively  simple  exterior  problem  using  the 
method  of  moments  to  determine  the  equivalent  excitation 
currents  and/or  fields  in  the  aperture  regions.  This  can  be 
done  independent  of  any  knowledge  of  the  interior  contents 
of  the  object  of  interest. 

Step  2  -  Use  of  the  FD-TD  method  to  analyze  the  relatively  complex 
interior  problem,  assuming  as  an  excitation  the  equivalent 
currents  or  fields  found  in  Step  1. 

In  this  way,  each  analysis  technique  is  applied  in  the  structure  size  and 
complexity  range  that  it  is  best  suited  for,  allowing  an  overall  solution 
that  is  accurate  for  large,  simple  structures  with  small,  complex  apertures 
or  other  fine  interior  details. 

IITRI  employed  the  hybrid  MOM/FD-TD  technique  to  solve  the  following 
electromagnetic  coupling  problems: 

Interacting  structure  -  Circular  (19.0  cm  diameter),  68.5  cm  long, 
open-ended  aluminum  cylinder  (shown  in  Figures  40,  41,  52,  and  53  of 
the  Final  Report); 

Incident  wave  -  300  MHz  plane  wave  propagating  towards  the  open  end 
of  the  cylinder.  Wave  angles  of  incidence  investigated: 

Case  1  -  Axial  incidence; 

Case  2  -  Oblique  incidence  (wavevector  at  a  45°  angle  with 

respect  to  the  cylinder  axis)  and  transverse  magnetic 
(TM)  polarization; 

Desired  fields  -  Each  component  of  total  E  and  total  H  in  the  axial 
cross  section  plane  of  the  cylinder  down  to  40  cm  from  the  open  end. 
First,  with  the  cross-section  plane  parallel  to  the  incident  H  component, 
and  again  with  the  plane  parallel  to  the  incident  E  component; 

Resolution  -  0.5  cm  uniformly  throughout  the  mapping  planes; 

Plotted  values  -  In  decibels  relative  to  an  incident  E"  of  1  volt/meter 


and  an  incident  H  of  1/377  ampere/meter. 


To  solve  these  coupling  problems,  IITRI  suitably  modified  two  existing 
FD-TD  computer  codes  to  accept  data  for  the  equivalent  aperture  electric  cur¬ 
rent  excitation.  These  data  were  then  obtained  from  Drs.  D.  Wilton  and 
A.  Glisson*  who  employed  their  MOM  code  for  the  electric  currents  over  the 
cylinder  aperture  locus  when  shorted,  as  specified  by  the  Schelkunoff  theorem. 
Details  of  the  models  are  given  in  Sections  7.4.1  and  7.4.2  of  the  Final 
Report. 

IITRI  also  employed  the  hybrid  MOM/FD-TD  technique  to  solve  the  electro¬ 
magnetic  coupling  problem  of  the  loaded  missile  guidance  section,  which  had 
been  previously  solved  by  IITRI  during  this  research  program  using  the  pure 
FD-TD  method  (Task  2,  Case  2).  The  excitation  consisted  of  MOM-derived  data 
for  the  equivalent  electric  current  excitation  over  the  nose  aperture  and  the 
sleeve-fitting  aperture  when  shorted,  as  specified  by  the  Schelkunoff  theorem. 
These  data  were  again  obtained  from  Drs.  D.  Wilton  and  A.  Glisson.  Details 
of  this  model  are  given  in  Section  7.4.3  of  the  Final  Report. 

All  hybrid  MOM/FD-TD  numerical  results  were  then  compared  to  available 
pure  FD-TD  results  (for  the  open-ended  cylinder  at  axial  incidence,  and  the 
loaded  missile  guidance  section),  and  to  available  pure  MOM  results  (for  the 
open-ended  cylinder  at  axial  incidence  and  at  oblique  incidence).  The  goal 
was  to  establish  the  consistency  and  accuracy  of  the  hybrid  modeling  tech¬ 
nique. 

4.0  RESULTS  AND  DISCUSSION 

Completion  of  Program  Task  1  demonstrates  the  capability  of  the  FD-TD 
method  to  accurately  model  the  electromagnetic  fields  coupled  into  a  cylin¬ 
drical  structure  which  is  illuminated  with  an  electromagnetic  plane  wave  hav¬ 
ing  an  arbitrary  polarization.  Elements  of  Task  5  demonstrate  this  modeling 
capability  for  an  arbitrary  wave  angle  of  incidence,  as  well. 

Completion  of  Program  Task  2  demonstrates  the  capability  of  the  FD-TD 
method  to  map  fields  coupled  into  a  complex  conical  structure  that  has  both 
irregular  shaped  apertures  and  interior  dielectric  and  metallic  materials. 

♦Consultants  to  IITRI  from  the  EE  Department,  University  of  Mississippi. 
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Determination  of  the  accuracy  bound  for  this  method,  as  applied  to  complex 
structures,  awaits  the  results  of  future  experimental  programs  since  other 
numerical  approaches  cannot  deal  directly  with  this  level  of  complexity.  No 
comparative  data  is  available  at  this  time. 

Completion  of  Program  Task  3  demonstrates  the  ability  of  the  FD-TD 
technique  to  directly  predict  the  current  that  the  coupled  fields  may  induce 
on  a  wire  or  bundle  of  wires  which  is  located  inside  a  complicated  structure 
(such  as  a  missile  guidance  section).  Here,  the  observed  parallelism  of  the 
computed  magnetic  field  contours  near  the  modeled  wire(s)  is  employed  in  com¬ 
bination  with  Ampere's  law  to  translate  FD-TD-computed  magnetic  fields  to 
wire  current  magnitudes.  Confirmation  of  the  FD-TD-computed  wire  currents  is 
required  in  a  future  experimental  program. 

Completion  of  Program  Task  4  demonstrates  the  computer  resources  needed 
to  implement  the  FD-TD  technique  for  a  given  size  structure  illuminated  by 
electromagnetic  waves  of  a  given  frequency.  This  task  leads  to  the  conclu¬ 
sions  that:  1)  If  on-line  mass  storage  is  used,  the  structure  size  is  limited 
mainly  by  the  maximum  "wall  clock  time"  of  the  program  which  is  acceptable  to 
both  the  user  of  the  program  and  the  operator  of  the  computing  center;  2)  If 
mass  storage  is  not  used,  the  structure  size  is  sharply  limited  by  the  avail¬ 
able  central  memory  because  larger  objects  (or  ones  with  higher  spatial  res¬ 
olution)  simply  cannot  be  specified  within  the  bounds  of  the  memory  field. 

Completion  of  Program  Task  5  leads  to  the  following  conclusions:  1) 
Hybrid  MOM/FD-TD  results  are  consistent  with  previous  data  derived  from  both 
pure  MOM  and  pure  FD-TD  methods,  implying  the  validity  of  the  hybrid  formu¬ 
lation;  2)  Arbitrary  angles  of  wave  incidence  can  be  treated  using  the 
hybrid  MOM/FD-TD  method  simply  by  specifying  the  proper  aperture  equivalent 
current  excitation,  without  having  to  change  the  FD-TD  geometry  of  the  struc¬ 
ture;  3)  Wires  or  other  structures  in  the  interior  of  a  cavity  strongly 
coupled  to  an  aperture  can  be  easily  and  consistently  modeled  using  the 
hybrid  MOM/FD-TD  approach;  and  4)  Complex  sub-sections  of  elongated  simply- 
shaped  structures  are  prime  candidates  for  detailed  modeling  of  the  interior 
penetrating  fields  via  the  hybrid  MOM/FD-TD  method. 
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Overall,  under  RADC  Contract  F30602-79-C-0039,  IITRI  demonstrated  that 
the  FD-TD  method  can  be  successfully  applied  to  electromagnetic  coupling 
problems  involving  highly-conducting  structures  with  hole  and  sleeve-type 
apertures  illuminated  by  a  plane  wave  having  an  arbitrary  polarization  and 
angle  of  incidence.  Further,  IITRI  showed  that  the  FD-TD  method  can  be  ap¬ 
plied  to  complex  cavity-like  structures  having  internal  metal  and  dielectric 
materials,  as  well  as  connecting  wires.  Accuracy  of  the  FD-TD  results  was 
very  good  relative  to  the  uncertainties  of  available  experimental  and  numer¬ 
ical-theory  approaches.  Convergence  of  the  electromagnetic  fields  to  the 
sinusoidal  steady  state  occurred  within  about  3  cycles  of  the  incident  wave 
when  a  slight  value  of  isotropic  conductivity  was  assigned  to  the  interior  of 
the  structures  modeled.  This  resulted  in  program  central  processor  times  of 
less  than  5  minutes  for  FD-TD  lattices  containing  as  many  as  1.5  million 
unknowns  time-stepped  to  n  ,  =  800  using  the  Control  Data  STAR-100  and 

Cyber  203  computers. 

During  this  research  program,  IITRI  also  established  the  feasibility  of 
a  hybrid  MOM  /  FD-TD  analysis  technique  based  upon  the  use  of  Schelkunoff ‘s 
equivalent  current  theorem  at  apertures  of  a  cavity.  This  hybrid  technique 
was  shown  to  give  consistent  results  for  structures  illuminated  at  both 
oblique  and  axial  incidence  by  plane  waves.  Further,  this  technique  was 
shown  to  give  consistent  results  for  apertures  strongly  coupled  to  internal 
wires  or  other  metal  and  dielectric  structures. 

IITRI  has  shown  that  the  pure  FD-TD  and  hybrid  MOM  /  FD-TD  methods  have 
great  promise  for  applications  involving  complex  conducting  and  dielectric 
structures  illuminated  by  plane  waves  at  arbitrary  angles  of  incidence  and 
polarization.  Further,  great  promise  is  shown  for  those  structures  that  are 
simultaneously  electrically  large  compared  to  a  wavelength  and  penetrated  by 
locally  complex  cavity-backed  apertures  having  dimensions  comparable  to  a 
wavelength.  Such  structures  requiring  resolution  of  both  large  and  small 
details  at  the  same  time  have  not  been  well  treated  by  any  one  previous 
analytical  or  numerical  approach. 
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No.  E6461 .  The  report  is  in  two  volumes.  Volume  1  covers  details  of  the 
technical  work,  including  relevant  theory  and  numerical  results.  Volume  2 
provides  listings  of  the  Fortran  computer  programs  used  to  obtain  the  results 
of  Volume  1 . 

The  principal  investigator  on  this  program  was  Dr.  Allen  Taflove,  with 
support  provided  by  Dr.  Korada  Umashankar  of  IITRI  on  Task  5;  Mr.  Neil 
Robertson  of  IITRI  on  data  reduction  and  contour  map  plotting;  and  Prof. 
Donald  Wilton  and  Dr.  Allen  Glisson  of  the  University  of  Mississippi  on  com¬ 
putation  of  frequency-domain  data  for  checking  of  results  and  for  Task  5. 

The  project  duration  was  7  December  1978  to  29  February  1980. 

Respectfully  submitted, 
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EVALUATION 


The  purpose  of  this  effort  was  to  further  examine  the  capabilities  of  the  Finite 
Difference-Time  Domain  (FD-TD)  technique  when  used  to  predict  electromagnetic 
coupling.  This  technique  has  demonstrated  capabilities  far  superior  to  any 
developed  thus  far.  The  FD-TD  technique  can  model  a  structure  that  consists  of  a 
number  of  materials  including  dielectrics  and  finite  conductors.  This  technique  has 
been  used  to  map  the  electromagnetic  fields  coupled  through  complex  shaped 
apertures  and  the  currents  induced  on  wires  behind  the  aperture.  There  is  no 
restriction  on  the  direction  of  incidence  or  polarization. 

The  FD-TD  technique  has  been  hybridized  with  "Method  of  Moments"  to 
improve  the  economy  of  computer  resources,  and  increase  the  range  of 
applicability  of  both  techniques. 

In  order  for  the  FD-TD  technique  to  receive  wide  distribution  it  is 
recommended  that  this  technique  be  developed  into  a  user  oriented  computer  code. 
As  FD-TD  is  used  more,  the  results  can  be  used  to  impact  further  development  of 
this  technique. 


DANIEL  E.  WARREN 
Project  Engineer 
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VOLUME  1  -  THEORY  AND  NUMERICAL  RESULTS 

1.0  INTRODUCTION 

Electromagnetic  penetration  problems  are  difficult  to  treat  with  many 
analytical  or  numerical  methods  because  of  the  inability  of  these  methods  to 
simply  deal  with  the  effects  of  structure  apertures,  curvatures,  corners,  and 
internal  contents.  Usually,  only  relatively  simple  geometries  and  apertures 
are  studied  in  an  attempt  to  gain  insight  into  the  key  penetration  mechanisms 
and  to  allow  an  indirect  estimate  of  the  penetration  for  mere  complicated 
problems. 

This  research  program  investigated  the  application  of  a  new  approach  for 
the  direct  modeling  of  electromagnetic  interaction  problems:  the  finite- 
difference,  time-domain  (FD-TD)  solution  of  Maxwell's  equations.  The  FD-TD 
method  treats  the  irradiation  of  a  structure  as  an  initial  value  problem.  At 
t  =  0,  a  plane  wave  source  of  frequency,  f,  is  assumed  to  be  turned  on.  The 
propagation  of  waves  from  this  source  is  simulated  by  solving  a  finite- 
difference  analog  of  the  time-dependent  Maxwell's  equations  on  a  lattice  of 
cells,  including  the  structure.  Time-stepping  is  continued  until  the  sin¬ 
usoidal  steady  state  is  achieved  at  each  cell.  The  field  envelope,  or  maxim¬ 
um  absolute  value,  during  the  final  half  wave-cycle  of  time-stepping  is  taken 
as  the  magnitude  of  the  phasor  of  the  steady-state  field. 

This  method  has  two  key  advantages  relative  to  available  modeling  ap¬ 
proaches.  First,  it  is  simple  to  implement  for  complicated  metal/dielectric 
structures  because  arbitrary  electrical  parameters  can  be  assigned  to  each 
lattice  cell  using  a  data  card  deck.  Second,  its  computer  memory  and  running 
time  requirement  is  not  prohibitive  for  many  complex  structures  of  interest. 
In  recent  work,  the  FD-TD  method  has  been  shown  capable  of  accurately  solving 
for  hundreds  of  thousands  of  unknown  field  components  within  a  few  minutes 
on  an  array-processing  computer.  Consistently,  a  +  1  dB  accuracy  relative  to 
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known  analytical  and  experimental  bench  marks  has  been  achieved  for  a  variety 
of  dielectric  and  metal  geometries. 

The  objective  of  the  present  study  is  to  more  completely  evaluate  the 
suitability  of  the  FD-TD  method  to  determine  the  amount  of  electromagnetic 
coupling  through  an  aperture  into  an  enclosed  conducting  container  and  the 
interaction  and  coupling  of  the  penetrating  fields  with  internal  electronics. 
Two  specific  container  models  are  used  for  the  evaluation.  The  first,  a  con¬ 
ducting  cylinder  with  one  open  end.  The  other,  the  guidance  section  of  a 
missile.  Each  of  these  two  configurations  is  modeled  to  calculate  the  elec¬ 
tromagnetic  field  coupled  into  the  structure. 

The  ultimate  aim  of  research  in  this  area  is  two-fold.  First,  develop 
an  easily  used  general  code  solving  for  the  fields  within  arbitrary  metal/ 
dielectric  structures  having  multiple  apertures,  while  assuring  a  known 
bounded  level  of  uncertainty.  And  second,  develop  a  more  sophisticated  in¬ 
tuitive  understanding  of  basic  wave  penetration  mechanisms  in  time  domain, 
such  as  transient  propagation  through  beyond-cutoff  cavity  interiors,  field 
build-up  at  edges,  and  convergence  to  the  sinusoidal  steady  state. 

This  report  first  reviews  in  Section  2.0  the  basic  theory  behind  the 
FD-TD  method,  referencing  extensively  technical  report  RADC-TR-78-142,  "Time 
Domain  Solutions  for  Electromagnetic  Coupling,"  the  final  technical  report 
on  Contract  F30602-77-C-0163  (the  forerunner  of  the  present  contract).  This 
section  also  briefly  summarizes  the  results  of  the  previous  contract. 

Sections  3.0  through  7.0  document  the  results  of  the  present  research 
effort  in  terms  of  the  following  five  tasks: 

Task  1:  Prediction  of  the  Coupling  into  an  Open-Ended  Cylinder; 

Task  2:  Prediction  of  the  Coupling  into  a  Missile  Guidance  Section f 

Task  3;  Prediction  of  the  Coupling  to  Wires  or  Cable  Bundles? 

Task  4;  Determination  of  the  Size  Limits  of  an  Object  that  can  be 

Modeled  by  the  FD-TD  Method;  and 

Task  5:  Investigation  of  Interfacing  the  FD-TD  Method  with  Other 
Analysis  Techniques. 
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Appendices  A  and  B  provide  theoretical  and  numerical  background  to  the  results 
of  Task  5.  These  appendices  discuss  the  basis  of  the  hybrid  technique  in¬ 
volving  the  Method  of  Moments  and  the  FD-TD  method. 

2 . 0  THE  FO-TD  METHOD 

This  section  will  review  the  followinq  topics  relevant  to  the  use  of  the 
FO-TD  method:  1)  the  basic  ideas  behind  the  method;  2)  computational 
details;  and  3)  past  usages  and  results  in  RADC  Contract  F30602-77-C-0163. 

2 .  1  Basic  Ideas 

2.1.1  Wave  Tracking 

The  FD-TD  method  is  a  direct  solution  of  Maxwell's  time-dependent  curl 
equations.  The  goal  is  to  model  the  propagation  of  an  electromagnetic  wave 
into  a  volume  of  space  containing  a  dielectric  or  conducting  structure.  By 
time  stepping,  i.e.,  repeatedly  implementing  a  finite-difference  analog  of 
the  curl  equations  at  each  cell  of  the  corresponding  space  lattice,  the  inci¬ 
dent  wave  is  tracked  as  it  first  propagates  to  the  structure  and  then  inter¬ 
acts  with  it  via  surface  current  excitation,  diffusion,  penetration,  and 
diffraction.  Wave  tracking  is  completed  when  the  desired  late-time  or  sin¬ 
usoidal  steady  state  behavior  is  observed  at  each  lattice  cell.  The  rationale 
for  this  procedure  is  that  it  achieves  simplification  by  analyzing  the  inter¬ 
action  of  the  wavefront  with  portions  of  the  structure  surface  at  a  given 
instant  in  time,  rather  than  attempting  a  simultaneous  solution  of  the  entire 
problem. 

Time  stepping  for  the  FD-TD  method  is  accomplished  by  an  explicit  finite- 
difference  procedure  due  to  Yee  [1].  For  a  cubic-cell  space  lattice,  this 
procedure  involves  positioning  the  components  of  F  and  H  about  a  unit  cell 
of  the  lattice  as  shown  in  Figure  1,  and  evaluating  E  and  FI  at  alternate  half 
time  steps.  In  this  manner,  centered  difference  expressions  can  be  used  for 
both  the  space  and  time  derivatives  to  attain  second-order  accuracy  in  the 
space  and  time  increments  without  requiring  simultaneous  equations  to  compute 
the  fields  at  the  latest  time  step. 
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Fig.  1  POSITIONS  OF  THE  FIELD  COMPONENTS 

ABOUT  A  UNIT  CELL  OF  THE  YEE  LATTICE 


4 


2.1.2  Modeling  the  Structure  of  Interest 

The  finite-difference  formulation  of  the  FD-TD  method  allows  the 
straightforward  modeling  of  the  surfaces  and  interiors  of  arbitrary  dielectric 
or  conducting  structures.  The  structure  of  interest  is  mapped  into  the  space 
lattice  by  first  choosing  the  space  increment  and  then  employing  a  data  card 
deck  to  assign  values  of  permittivity  and  conductivity  to  each  component  of  E. 
No  special  handling  of  electromagnetic  boundary  conditions  at  media  interfaces 
is  required  because  the  curl  equations  generate  these  conditions  in  a  natural 
way  by  themselves.  Therefore,  the  basic  computer  program  need  not  be  modified 
to  change  from  structure  to  structure.  In  this  manner,  inhomogeneities  or 
fine  details  of  the  structure  can  be  modeled  with  a  maximum  resolution  of  one 
unit  cell;  thin  surfaces  can  be  modeled  as  infinitely  thin,  stepped-edge 
sheets. 

The  explicit  formulation  of  the  FD-TD  method  is  particularly  suited  for 
programming  with  minimum  storage  and  execution  time  using  recently  developed 
array-processing  computers.  First,  the  required  computer  storage  and  running 
time  increases  linearly  with  N,  the  total  number  of  unknown  field  components. 
Computer  techniques  (such  as  the  method  of  moments)  which  require  the  solution 
of  simultaneous  equations  usually  have  a  storage  requirement  proportional  to 
N2  and  a  running  time  proportional  to  N2  to  N3  [2].  Second,  since  all  FD-TD 
operations  are  explicit  and  can  be  performed  in  parallel,  rapid  array-proces¬ 
sing  techniques  can  be  readily  applied.  As  will  be  demonstrated  later,  these 
can  be  employed  to  solve  for  10s  -  106  field  components  in  a  single  FD-TD 
problem,  as  opposed  to  a  maximum  of  about  103  field  components  for  convention¬ 
al  approaches  using  simultaneous  equation  solutions  [2]. 

2.1.3  Algorithm  Refinements 

Yee  applied  the  FD-TD  method  to  compute  the  waveforms  of  TM  and  TE  pulses 
scattered  from  infinitely  long,  rectangular  cross-section,  conducting  cylin¬ 
ders.  Other  workers  investigated  electromagnetic  pulse  interactions  in  time- 
varying  inhomogeneous  media  [3],  with  metallic  bodies  of  revolution  [4],  and 
with  detailed  models  of  aircraft  [5].  Three  distinct  problems  emerged  in  the 
process  of  adapting  the  FD-TD  method  to  model  realistic  situations: 
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1.  Lattice  truncation  conditions.  The  field  components  at  the  lattice 
truncation  planes  cannot  be  determined  directly  from  the  Maxwell's 
equations  analog  and  must  be  computed  using  an  auxiliary  radiative 
truncation  condition.  However,  care  must  be  exercised  because  this 
condition  must  not  cause  excessive  spurious  reflection  of  waves  scat¬ 
tered  outward  by  the  structure  modeled.  The  goal  is  to  formulate 
truncation  planes  as  close  as  possible  to  the  structure  (to  minimize 
computer  storage),  and  yet  achieve  virtual  invisibility  of  these  planes 
to  all  possible  waves  within  the  lattice. 

2.  Plane  wave  source  condition.  The  simulation  of  either  an  incident 
plane-wave  pulse  or  single-frequency  plane  wave  should  not  take  exces¬ 
sive  storage  or  cause  spurious  wave  reflections.  The  former  would 
occur  if  the  incident  wave  is  programmed  as  an  initial  condition;  the 
latter  would  occur  if  the  incident  wave  is  programmed  as  a  fixed  field 
excitation  along  a  single  lattice  plane. 

3.  Sinusoidal  steady  state  information.  Such  data  can  be  obtained  either 
by  (1)  directly  programming  a  single-frequency  incident  plane  wave  or 
(2)  performing  a  separate  Fourier  transformation  step  on  the  pulse  wave¬ 
form  response.  Both  methods  require  time  stepping  to  a  maximum  time 
equal  to  several  wave  periods  at  the  desired  frequency.  The  second 
method  has  two  additional  requirements.  First,  a  short  rise  time  pulse 
suffers  from  accumulating  waveform  error  due  to  overshoot  and  ringing  as 
it  propagates  through  the  space  lattice.  This  leads  to  a  numerical  noise 
component  which  should  be  filtered  before  Fourier  transformation. 

Second,  Fourier  transformation  of  many  lattice-cell  field  vs.  time  wave¬ 
forms  (each  probably  extending  over  hundreds  of  time  steps)  would  sig¬ 
nificantly  add  to  the  total  requirements  for  computer  storage  and  execu¬ 
tion  time. 

Previous  work  leading  to  RADC  Contract  F30602-77-C-01 63  described  efforts 
to  solve  the  above  problems  for  the  case  of  the  cubic-unit-cell  space  lattice 
[6],  [7].  Simple  truncation  conditions  were  developed  for  two-  and  three- 
dimensional  lattices  that  reduced  the  reflection  coefficient  of  closely- 
positioned  truncation  planes  to  the  order  of  0.1  for  waves  of  arbitrary  incid¬ 
ence.  A  plane-wave  source  condition  was  described  that  allowed  generation  of 
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an  arbitrary  pulsed  or  sinusoidal  incident  wave  without  requiring  any  addi¬ 
tional  storage  and  without  causing  spurious  wave  reflections.  Finally,  it  was 
shown  that  sinusoidal  steady-state  data  could  be  computed  using  the  FD-TD 
method  by  directly  programming  a  single-frequency  incident  plane  wave  and 
time-stepping  to  the  steady  state  over  several  periods  of  the  incident  wave. 
Observed  accuracy  for  low-permittivity  cylindrical  dielectric  geometries  was 
±  1  dB  at  points  of  electric  field  maxima,  and  +  1  lattice  cell  in  locating 
electric  field  maxima  and  minima. 


Under  RADC  Contract  F30602-77-C-0163 ,  additional  FD-TD  method  provisions 
were  developed  to  allow  accurate  solutions  for  the  fields  penetrating  metal 
structures  having  one  or  more  apertures  [8].  This  class  of  structures  ex¬ 
hibits  a  different  electromagnetic  penetration  response  than  convex  dielectric 
objects  due  to  such  phenomena  as  aperture  coupling,  edge  effects,  and  beyond- 
cutoff  interiors.  Further,  the  reflection  coefficient  of  the  metal  surfaces 
modeled  is  larger  than  that  of  dielectric  surfaces,  leading  to  increased  dif¬ 
fraction  and  the  necessity  for  improved  lattice  truncation  conditions.  The 
new  FD-TD  method  provisions  are  now  qualitatively  discussed. 


The  first  algorithm  development  involves  the  introduction  of  an  aniso¬ 
tropic  lossy  air  medium  outside  of  the  modeled  structure  to  aid  in  reducing 
the  effective  reflection  coefficient  of  the  lattice  truncation  planes. 
Modeling  of  such  a  medium  is  possible  because  separate  difference  equations 
are  available  for  each  field  vector  component,  thus  allowing  individual  spec¬ 
ification  of  loss  to  each  field  component.  For  the  three-dimensional  case, 
it  is  useful  to  specify  an  anisotropic  electric  field  loss,  ogxt,  in  the 
free-space  region  exterior  to  the  structure.  This  loss  term  is  used  only  in 
the  difference  equations  for  the  electric  field  components  not  present  in  the 
incident  plane  wave.  It  is  selected  to  induce  1/e  decay  of  these  components 
as  they  propagate  over  a  distance  about  equal  to  the  total  dimension  of  the 
lattice.  These  choices  avoid  any  exponential  decay  of  the  incident  wave, 
which  would  lead  to  serious  modeling  error  for  structure  sizes  comparable  to 
the  1/e  distance.  Further,  the  resulting  modified  scattered  wave  is  reduced 
in  amplitude  upon  reaching  the  lattice  truncations  and  is  further  damped  as 
its  remnants  reverberate  between  the  structure  and  the  truncations.  Applica¬ 
tion  of  this  development  has  been  found  to  be  important  in  maintaining  a  ±  1 


dB  accuracy  of  the  FD-TD  results  for  metal  structures.  For  low  permittivity 
dielectric  structures,  the  assumption  of  completely  lossless  air  in  the  exter¬ 
ior  region  is  sufficient  for  this  accuracy  level. 


The  second  algorithm  development  involves  the  introduction  of  a  small 
isotropic  conductivity,  o.^,  within  the  interior  of  the  modeled  metal  struc¬ 
ture.  The  purpose  of  modeling  the  slightly  lossy  air  is  to  cause  the  rever¬ 
berating  fields  within  the  structure  to  converge  more  rapidly  to  the  expected 
steady-state  condition,  especially  for  the  case  when  the  interior  cavity  vol¬ 
ume  is  normally  beyond  cutoff  at  the  frequency  of  the  incident  plane  wave. 

In  this  case,  the  transient  internal  fields  are  highly  reactive  (carry  little 
real  power  flow)  and  quickly  dissipate  when  forced  to  supply  energy  to  main¬ 
tain  an  electric  field  distribution  across  a  slightly  lossy  medium. 

Accurate  results  are  obtained  in  these  circumstances  when  cr.  .  is 

l  nt 

chosen  to  provide  an  equivalent  10  dB  of  loss  for  a  plane  wave  propagating 
the  full  length  of  the  metal  structure  interior. 

2.2  Computational  Details  for  a  Uniform,  Cubic  Lattice 

2.2.1  System  of  Finite-Difference  Equations 

Using  the  MKS  system  of  units,  and  assuming  that  the  dielectric  param¬ 
eters,  y,  e,  and  a,  are  independent  of  time,  the  following  system  of  scalar 
equations  is  equivalent  to  Maxwell's  equations  in  the  rectangular  coordinate 
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Yee  [1]  originally  introducteci  a  set  of  finite-difference  equations  for  the 
system  of  Equations  (la)  -  (If).  Following  Yee's  notation,  we  denote  a  space 
point  in  a  cubic  lattice  as 

(i  ,j  ,k)  =  (i6,j6,k6)  (2) 

and  any  function  of  space  and  time  as 

Fn(i,j,k)  =  F(i6,j<5,k6,n6t) ,  (3) 

where  5  =  6x  =  6y  =  6z  is  the  space  increment,  6t  is  the  time  increment,  and 
i,  j,  k,  and  n  are  integers.  Yee  used  centered  finite-difference  expressions 
for  the  space  and  time  derivatives  that  are  both  simply  programmed  and 
second-order  accurate  in  6  and  in  <5t,  respectively: 
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To  achieve  the  accuracy  of  Equation  (4),  and  to  realize  all  of  the  space 
derivatives  of  Equations  (la)  -  (If),  Yee  positioned  the  components  of  E  and 
H  about  a  unit  cell  of  the  lattice  as  shown  in  Figure  1.  To  achieve  the  ac¬ 
curacy  of  Equation  (5),  he  evaluated  E  and  F  at  alternate  half  time  steps. 
The  result  of  these  assumptions  is  the  following  general  system  of  finite- 
difference  equations  for  the  system  of  Equations  (la)  -  (If): 
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With  the  system  of  Equations  (6a)  -  (6f),  the  new  value  of  a  field  vec¬ 
tor  component  at  any  lattice  point  depends  only  on  its  previous  value  and  on 
the  previous  values  of  the  components  of  the  other  field  vector  at  adjacent 
points.  Therefore,  at  any  given  time  step,  the  computation  of  a  field  vector 
may  proceed  one  point  at  a  time. 

Many  electromagnetic  interaction  problems  involve  nonpermeable  media  and 
can  be  approached  using  a  fixed  time  step  and  space  increment.  For  such  prob¬ 
lems  (including  the  cylinder  and  nose  cone  geometries  specified  for  this  re¬ 
search  effort),  the  quantity  6t/v(i,j,k)6  is  constant  for  all  (i,j,k)  of  the 
lattice,  and  the  Yee  system  of  Equations  (6a)  -  (6f)  can  be  simplified  to 
reduce  computer  running  time  in  the  following  manner.  We  define  the  con¬ 
stants  : 
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R  »  6t/2e0 

(7a) 

Ra»  6t2/(62w0e0) 

(7b) 

Rb  6t/wo6 

(7c) 

1  -  Ro(m)/er(m) 

(7d) 

Ca^  1  +  Ro(m)/er{mF 

Ra 

Cb^  "  er(m)  +  RolmJ 

(7e) 

where  m  is  an  integer  denoting  a  particular  dielectric  or  conducting  medium 
in  the  space  to  be  modeled.  We  also  define  the  proportional  electric-field 
vector 

E  -  Rb  F  (8) 


Using  the  definitions  of  Equations  (7a)  -  (7e)  and  (8),  we  rewrite  Equations 
(6a)  -  (6c)  as: 


n+% 

Hx(i ,j+%,k+%) 


n-%  „n 

Hx(i,j+%,k+%)  +  Ey ( i »j+%,k+l )  -  Ey(i ,j+%,k)  + 
E"(i.j.k+%)  -  E*(i ,j+l,k+%) 


(9a) 


n+i; 

Hy(i+%.j,k+%) 


n-% 

Hy(i+%.j,k+%) 
Ex(i+%,j,k)  - 


+  Ez(i+U»k+%)  -  Ez(i,j,k+%)  + 
Ex(i+%,j,k+l) 


(9b) 


n+% 

Hz(i+%J+%,k) 


Hz(t+%,j+%,k)  +  Ex(i+%,j+l.k)  -  Ex(i+%,j,k)  + 


Ey(1 »j+%»k)  -  Ey(i+1 ,j+%,k) 


(9c) 
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This  modification  eliminates  the  three  multiplications  needed  by  Yee  in  the 
H  part  of  the  algorithm.  Further,  we  rewrite  Equations  (6d)  -  (6f)  as: 


m  =  MEDIA(i+%,j,k) 

n+1  ~n  n+^ 

Ex(i+%,j,k)  =  Ca(m)Ex(i+%,j,k)  +  Cb(m)[Hz(i+%,j+%,k)  - 

n+%  o+% 

.k-%)  -  Hy(i+%,j ,k+%)] 


m  =  MEDIA(i,j+%,k) 


n+1  n  n+% 

E,(i,j+%,k)  =  C  (m)E  (i ,j+%,k)  +  Cb(m)[H  (i ,j+%,k+%) 

y  ay  u  a 

n+%  n+% 

Hz(i-%,j+%,k)  -  Hz(i+%,j+%,k)] 


m  =  MEDI A( i ,j ,k+%) 


n+1  n  n+% 

Ez(i,j,k+%)  =  Ca(m)Ez(i,j,k+%)  +  Cb(m)[Hy(i+%,j,k+%) 


n+%  n+% 

Hx(i ,j-%,k+%)  -  Hx ( i , j+%, k+%) ] 


n+% 

Hz(i+%»j-%.k)  + 
(9d) 


n+% 

Hx(i ,j+%,k-%)  + 
(9e) 


n+% 

-  Hy(i-%,j ,k+%)  + 
(9f) 
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This  modification  eliminates  the  need  for  computer  storage  of  separate  t  and 
a  arrays.  Now,  only  a  MEDIA  array  which  specifies  the  type-integer  of  the 
dielectric  or  conducting  medium  at  the  location  of  each  electric  field  compon¬ 
ent  in  the  lattice  need  be  stored.  In  addition,  the  c  and  o  of  each  medium 

r* 

can  now  be  changed  without  having  to  re-punch  a  large  data  card  deck,  if  the 
basic  structure  geometry  is  unchanged.  Such  a  change  involves  only  the  re¬ 
calculation  of  the  few  values  of  C,(m)  and  C.(m). 

d  D 

2.2.2  Choice  of  Space  and  Time  Increments 

The  choice  of  5  and  fit  is  motivated  by  the  reasons  of  accuracy  and 
algorithm  stability,  respectively.  To  insure  the  accuracy  of  the  computed 
spatial  derivatives  of  the  electromagnetic  fields,  fi  must  be  small  compared 
to  a  wavelength  (usually  <_  A/10).  Further,  to  insure  that  the  cubic  lattice 
approximation  to  the  surfaces  of  the  structure  modeled  is  not  too  coarse,  fi 
must  be  small  compared  to  the  overall  dimensions  of  the  structure. 


To  insure  the  stability  of  the  time-stepping  algorithm  of  Equations 
(9a)  -  (9f),  fit  is  chosen  to  satisfy  the  inequality 


fit 


( 


fix' 


_L 

6y2 


-T> 

fiz 


-1 

cmax 


(10) 


<  _ §. -  (for  a  cubic  lattice) 

c  /~3 
max 

where  c,  is  the  maximum  wave  phase  velocity  within  the  model.  The  corres- 
max 

ponding  stability  criterion  set  forth  by  Yee  in  Equations  (7)  and  (8)  of  his 
paper  is  incorrect.  The  derivation  of  Equation  (10)  was  outlined  in  Section 

3.4.2  of  Reference  8. 


2.2.3  Lattice  Truncation  Conditions 

A  basic  consideration  with  the  FD-TD  lattice  is  the  treatment  of  the 
field  vector  components  at  the  lattice  truncation  planes.  Inspection  of 
Equations  (9a)  -  (9f)  indicates  that  the  values  of  such  components  cannot  be 
determined  from  the  system  of  finite-difference  equations  because  of  the 
centered  nature  of  the  spatial  derivatives.  Therefore,  these  values  must  be 
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computed  using  an  auxiliary  truncation  condition.  However,  great  care  must  be 
taken  because  this  condition  must  not  cause  the  spurious  reflection  of  waves 
scattered  outward  from  the  structure  modeled,  as  observed  by  Yee.  The  goal 
of  formulating  the  truncation  condition  is  to  make  the  lattice  truncation 
planes  invisible  to  all  possible  waves  propagating  within  the  lattice,  as 
shown  in  Figure  2. 

A  desirable  truncation  condition  relates  in  a  simple  way  the  values  of 
the  field  components  at  the  truncation  planes  to  field  component  values  at 
points  one  or  more  6  within  the  lattice.  We  now  consider  examples  of  such  a 
truncation  for  cases  of  FD-TD  lattices  in  one  and  three  dimensions. 

One-Dimensional  Case 

For  simplicity,  we  consider  waves  having  only  the  Ez  and  components 
and  propagating  in  the  ±y  directions.  The  one-dimensional  FD-TD  lattice  is 
simply  a  y-directed  line  of  points  having  the  Ez  and  H  components  inter¬ 
leaved  and  separated  from  each  other  by  0.5  6y.  The  lattice  is  assumed  to 
extend  from  an  Ez  component  at  point  y  =  0  to  another  Ez  component  at  point 
y  =  Jfiy.  A  time  step  of  fit  =  5y/c  is  used;  a  value  which  is  the  maximum 
allowed  by  the  stability  condition  of  Equation  (10)  for  this  lattice  (fix  = 
fiz  =  °°). 

Subject  to  these  assumptions,  the  truncation  condition  at  point  y  =  0, 

Ez  (0)  =  (1)  (11a) 

simulates  the  free  space  propagation  of  the  magnitude  of  Ez  from  the  point 
"1"  to  the  truncation  point  "0"  in  one  time  step  (the  free-space  propagation 
delay  implied  by  the  time-step  relation).  This  is  an  exact  truncation  for 
this  lattice  in  that  all  possible  -y-directed  waves  are  absorbed  at  0  without 
reflection.  If  we  wish  to  simulate  the  truncation  of  the  lattice  at  point 
y  =  Jfi,  the  truncation  condition 

e"  (J)  =  e"_1  (0-1)  (11b) 

is  exact  for  all  possible  +y-directed  waves  at  this  point. 
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Lattice  truncation  plane 


Fig.  2  IDEAL  FD-TD  LATTICE  TRUNCATION  CONDITIONS 


Three-Dimensional  Case 


Here,  we  consider  waves  having  all  six  field  components  and  propagating 
in  all  possible  directions.  The  lattice  is  assumed  to  extend  from: 

Hy  and  Hy  components  at  x  =  ^6  to 

H  and  components  at  x  =  (I  +  '-2)6; 

Ex  and  components  at  y  =  0  to 

and  E?  components  at  y  =  J6; 

E^  and  E^  components  at  z  -  0  to 

Ex  and  E  components  at  z  -  K6. 

A  time  step  of  <St  =  6/2 c  is  used,  a  value  which  is  about  13%  lower  than  the 
maximum  allowed  (6t  =  5 ! c)  by  the  stability  condition  of  Equation  (10)  for 
this  lattice  (6x  =  6y  =  5z  =  6). 

No  simple,  exact  truncation  condition,  analogous  to  Equations  (11a)  and 
(lib),  is  apparent  for  this  three-dimensional  space  lattice.  This  is  because 
we  cannot  assume  the  outgoing  waves  to  be  plane  and  normally  incident  on  one 
lattice  boundary.  At  any  truncation  point,  the  local  angle  of  incidence  of 
these  waves  relative  to  the  truncation  plane  is  unknown.  Further,  several 
different  waves  having  different  local  angles  of  incidence  may  arrive  at  the 
same  time.  No  simple  truncation  condition  can  account  for  all  of  these 
possibilities.  Therefore,  we  can  arrive  at  only  an  approximate  condiction 
that  reduces  the  effective  lattice  boundary  reflection  coefficient  to  an 
acceptable  level. 

A  set  of  simple,  approximate  truncation  conditions  that  has  been  used 
with  good  results  [8]  is  as  follows. 
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At  x  =  %  6: 
n 

Hy(%,j,k+%) 

n 

H2(%,j+%,k) 


n-2  n-2 

[Hy  (3/2 ,j ,k-%)  +  Hy  (3/2 ,j ,k+%) 


(12a) 


n-2 

+  H  (3/2 ,j ,k+3/2)]/3 


n-2  n-2 

[H2  (3/2,j+%,k-l)  +  Hz  (3/2 ,j+%,k) 


n-2 


(12b) 


+  Hz  (3/2 , j+%,k+l )]/3 


At  x  =  (I  +  %)  6: 
n 

Hy(I+%.j.k+%)  = 
n 

Hz(I+%.J+%»k)  = 


,  (12c) 

n-2  n-2  n-2 

[H  (I-%,j,k-%)  +  Hy  (I-%,j,k+%)  +  Hy  (I-%J,k+3/2)]/3 

(1 2d ) 

n-2  n-2  n-2 

[Hz  (I-%,j+%.k-l)  +  Hz  (I-%,j+%,k)  +  Hz  (I-%»j+%»k+l )]/3 


At  y  =  0: 

.n 

Ex(i+%,0,k) 

„n 

E2(i,0,k+%) 

At  y  =  J  6: 
n 

Ex(i+%,J,k) 

E2(i,J,k+%) 

At  z  =  0: 

.n 

ExO+%,j,0) 

yu+%,0) 


.n-2 

Ex  (i+%,1 *k) 
.n-2 

Ez  (i,l,k+%) 


(13a) 

(13b) 


n-2 

Ex  (i+%,J-l »k)  (13c) 


.n-2 

E2  (i ,J-1 ,k+%) 


(13d) 


_n-2 


.n-2 


(i-%.J,l)  +  Ex  (i+%,j,l )  +  e"  (i+3/2,j , 


>1 )]/3 


_~n-2  .n-2  n-2 

LE  (i-l.j+%,1)  +  E  (i,j+%,l)  +  E  (i+1  ,j+%,l  )]/3 


04a) 

(14b) 
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At  z  =  K  6: 

„n  n-2  n-2  n-2  04c) 

fx(i+%,j,K)  =  [Ex  ,K-1 )  +  Ex  (i+%,j,K-l)  +  Ex  (i+3/2 ,j ,K-1 )]/3 

-n  ^n-2  ^n-2  _n-2  (14d) 

Ey(i»j+%.K)  =  [Ey  ( i-1 ,j+%,K-l )  +  Ey  (i,j+%,K-l)  +  Ey  (i+1 ,j+%,K-l )]/3 

Equations  (12)  -  (14)  allow  the  field  value  at  any  truncation  point  to 
rise  to  approach  the  field  value  of  any  outgoing  wave,  thus  lowering  the  ef¬ 
fective  truncation  plane  reflection  coefficient.  This  is  done  by  modeling  the 
propagation  of  an  outgoing  wave  from  the  lattice  plane  adjacent  to  the  trunca¬ 
tion,  in  two  time  steps  (the  free-space  propagation  delay  implied  by  the  time- 
step  relation).  The  averaging  process  is  used  to  take  into  account  the  pos¬ 
sible  local  angles  of  incidence  of  the  outgoing  wave  at  the  truncation  and 
possible  multiple  incidences. 

Truncation  conditions  (12)  -  (14)  are  useful  for  an  assumed  +y-directed 
incident  plane  wave  with  field  components  E z  and  Hx>  For  such  a  wave.  Equa¬ 
tion  (13)  represents  exact  truncations  similar  to  Equation  (11).  In  addition. 
Equations  (12)  and  (14)  have  no  effect  on  the  propagation  of  such  a  wave, 

which  lacks  H  ,  H  ,  E  ,  and  E  .  Thus,  this  set  of  truncation  conditions  ef- 
y  z  x  y 

fectively  makes  the  lattice  boundary  planes  invisible  to  a  +y-directed  incid¬ 
ent  plane  wave. 


Use  of  Exterior-Region  Anisotropy 

One  way  of  reducing  spurious  reflections  at  the  lattice  truncations  is 
to  introduce  an  anisotropic  lossy  medium  outside  of  the  modeled  structure. 
Properly  constituted,  the  medium  would  attenuate  field  components  present  only 
in  the  scattered  wave,  leaving  the  incident  plane  wave  unaffected.  For  the 
three-dimensional  case,  this  can  be  easily  done  by  specifying  an  anisotropic 
conductivity,  aex^. ,  in  the  free-space  region  exterior  to  the  structure.  Equa¬ 
tion  (9f),  the  finite-difference  equation  for  Ez,  requires  o  =  0  to  insure 
that  the  incident  wave  is  not  attenuated.  However,  we  may  assume  a  small 
value  of  agxt  for  Equations  (9d)  and  (9e),  the  finite-difference  equations  for 

E  and  E  ,  without  affecting  the  propagation  of  the  incident  wave  or  the  pen- 
x  y 

etrating  wave  within  the  structure.  This  assumption  results  in  attenuation  of 

the  E  and  E  components  of  the  exterior  diffracted  wave,  and  thus,  reduces 
x  y 

z-directed  wave  reflections  at  the  lattice  truncations. 
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Effect  on  Algorithm  Stability 


The  stability  condition  of  Equation  (10)  is  valid  for  the  Yee,  or  null 
choice,  of  lattice  truncation  conditions.  This  is  because  Yee's  set  of  trun¬ 
cation  conditions  causes  total  reflection  of  all  lattice  wave  modes  at  the 
surface  planes  of  the  lattice,  and  thus,  introduces  no  new  wave  modes.  How¬ 
ever,  introduction  of  Equations  (12)  -  (14)  to  the  three-dimensional  algorithm 
is  found  to  increase  the  strictness  of  the  stability  condition.  For  the 
three-dimensional  case,  some  care  must  be  taken  to  avoid  algorithm  instabil¬ 
ity. 

The  nature  of  the  instability  of  the  three-dimensional  algorithm  is  of 
importance.  First,  it  is  late  in  appearance,  requiring  more  than  five- 

4 

hundred  time  steps  for  a  2  x  10  -  cell  lattice,  and  more  than  one-thousand 

4 

time  steps  for  a  6  x  10  -  cell  latice.  Second,  its  initial  visibility  is 

delayed  by  either  increasing  the  size  of  the  lattice,  or  by  increasing  the 
losses  of  the  dielectric  media  of  the  lattice.  This  suggests  the  importance 
of  wave  propagation  effects  in  the  growth  of  the  instability. 

There  are  two  likely  solutions  to  the  problem  of  algorithm  instability. 
First,  6t  can  be  reduced.  This,  however,  would  complicate  the  programming  of 
the  truncation  conditions  because  a  wave  would  no  longer  propagate  across  a 
free-space  unit  cell  of  the  lattice  in  an  integral  number  of  time  steps. 

The  second  solution  is  much  simpler  since  it  does  not  require  programming 
for  interpolating  field  values  at  the  truncation  planes  between  time  steps. 

5 

This  solution  is  merely  to  set  a  lower  bound  of  about  10  cells  on  the  size  of 
the  FD-TD  lattice  used,  thus  delaying  the  onset  of  the  instability  for  several 
thousand  time  steps.  In  most  cases,  such  a  delay  would  be  sufficient  for  the 
computed  solution  to  reach  the  sinusoidal  steady  state.  For  many  problems, 

5 

the  use  of  10  -  cell  lattice  is  not  at  all  extravagant  and  allows  the  problem 
of  algorithm  instability  to  be  essentially  forgotten.  The  computer  runs  in 
RADC  Contract  F30602- 77- C-01 63  for  the  cylinder  and  the  missile  nose  cone 
employed  this  solution.  They  used,  respectively,  94,000  cells  --  800  time 
steps  and  58,000  cells  --  900  time  steps,  without  any  apparent  instability  of 
the  computed  solution. 
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2.2.4  Plane  Wave  Source  Condition 


Another  basic  consideration  with  the  FD-TD  method  is  the  simulation  of 
the  continuous,  sinusoidal,  incident  plane  wave.  Yee  specified  the  shape  and 
direction  of  propagation  of  an  incident  wave  pulse  by  inserting  all  of  its 
field  values  as  initial  conditions  over  a  portion  of  the  lattice.  However, 
the  Yee  approach  is  clearly  inadequate  for  a  continuous  wave  train  because  a 
very  elongated  lattice  would  be  needed  to  contain  the  wave  as  an  initial  con¬ 
dition,  wasting  much  computer  storage. 

In  this  subsection,  we  discuss  the  simulation  of  an  incident,  +y-direct- 
ed  plane  wave  using  a  source  condition  localized  at  only  one  lattice  plane, 
and  invisible  to  all  scattered  waves  propagating  within  the  lattice.  This 
allows  a  compact  lattice  and  maximum  utilization  of  the  available  computer 
storage. 

The  most  simple  approach  to  this  problem  is  to  vary  the  electric  field 
at  all  points  along  the  lattice  plane  y  =  0  in  a  sinusoidal  manner.  This 
plane  would  then  radiate  the  desired  plane  wave.  However,  such  a  specifica¬ 
tion  of  field  values  at  a  lattice  truncation  plane,  without  consideration  of 
the  values  of  the  fields  of  any  possible  outgoing,  scattered  waves,  would 
cause  undesired  wave  reflections. 

A  more  desirable  plane  wave  source  condition  would  take  into  account  the 
scattered  fields  at  the  source  plane.  For  the  three-dimensional  case,  a  use¬ 
ful  wave  source  condition  at  plane  y  =  j$5  (near  y  =  0)  is  as  follows: 

-.n  ~n 

Ez(i,js^)  Rb  s i n ( 2rrf n6t )  +  Ez(i  ,j$ Mh)  (15) 

where  f  is  the  irradiation  frequency  and  is  defined  by  Equation  (7c). 
Equation  (15)  is  a  modification  of  the  Maxwell's  equations  algorithm  for  all 
points  on  the  lattice  plane  y  =  j$6.  At  each  point  on  this  source  plane,  the 
computer  first  calculates  E^  in  the  normal  manner  of  the  alogrithm,  and  stores 
the  value  in  memory.  Then,  the  value  of  the  sinusoid  is  calculated  and  added 
to  the  stored  value  of  e".  Finally,  this  modified  value  of  E^  is  stored  in 
memory.  In  effect.  Equation  (15)  simulates  the  linear  superposition  of  a 
+y-directed  plane  wave  and  the  ambient  field  along  the  source  plane.  This 


condition  permits  any  scattered,  outgoing  wave  to  propagate  right  through  the 
wave  source  plane  without  reflection,  and  reach  the  lattice  truncation  at 
y  =  0  to  be  absorbed. 

2.2.5  Symmetry  Conditions 

An  important  savings  of  computer  memory  and  program  execution  time  re¬ 
sults  if  even  symmetry  of  the  modeled  structure  about  one  or  two  lattice 
planes  can  be  assumed.  In  this  subsection,  we  discuss  the  programming  of 
this  symmetry  for  the  three-dimensional  case. 

For  the  three-dimensional  case,  the  modeled  structure  is  assumed  to  be 
evenly  symmetric  about  lattice  planes  x  =  (1+^)6  and  z  =  KS: 

e,cr(I+%+h,j  ,k)  =  £,o(T+%-h, j ,k) ,  U  =  U0  (16a) 

e,a(i,j,K+h)  =  e,a(i,j,K-h)  (16b) 

The  incident  radiation  is  assumed  to  be  a  +y-directed  plane  wave,  with  the 
field  components  Ez  and  Hx  naturally  having  even  symmetry  about  any  lattice 
plane  x  =  constant  or  z  =  constant.  Therefore,  we  conclude  that  the  E^  and 
Hx  components  of  the  total  field  possess  even  symmetry  about  the  lattice 
planes  x  *  (T+%) 6  and  z  =  (OS: 

~n  _  n 

E2(I+%+h,j,k+%)  =  Ez(I+%-h,j,k+%)  (17a) 


n+%_  n+%_ 

Hx  ( I+%+h,j+%,k+%)  =  Hx  (I+%-h,j+%,k+%)  (17b) 


n  _  n  _ 

Ez(i,j,K+h)  =  Ez(i,j,K-h) 


(17c) 


n+% 

Hx  (i,j+%,K+h) 


n+% 

Hx  (iJ+i,K-h) 


(17d) 
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To  develop  a  convenient  set  of  symmetry  conditions,  we  follow  a  procedure 
that  has  been  detailed  in  Section  3.4.5  of  Reference  8.  This  results  in 

n+1  n+%_  n+%_ 

Ex  (I+%,j,k)  «  Hy  (I+%,j,k+%)  =  Hz  (I+%,j+%,k)  =  0  for  all  n;  (18a) 

n+%  n  _  _n  _ 

Hz  (i+%,j+%,K)  =  Ex(i+%,j,K)  =  Ey(i,j+%,K)  =  0  for  all  n.  (18b) 

Equations  (18a)  and  (18b)  are  sufficient  to  truncate  the  FD-TD  lattice  at 
planes  x  =  (I+4)<5  and  z  =  K6,  respectively,  by  permitting  the  calculation  of 
the  complete  set  of  field  components,  with  full  specification  of  the  assumed 
even  symmetry. 

2.3  Initial  Evaluation  of  the  FD-TD  Method  - 

Results  of  Previous  RADC  Contract  F30602-77-C-0163 

This  section  first  summarizes  the  results  of  RADC  Contract  F30602-77-C- 
0163  previously  carried  out  to  evaluate  the  suitability  of  the  FD-TD  method 
to  determine  electromagnetic  coupling  through  an  aperture  into  an  enclosed 
conducting  container.  This  is  then  followed  by  a  discussion  of  the  basis  of 
the  present  research  to  more  fully  evaluate  the  FD-TD  method. 

2.3.1  Objective 

The  FD-TD  method  allows,  in  principle,  the  computation  of  the  internal 
EM  fields  of  complex  conducting  geometries.  However,  prior  to  the  work  on 
RADC  Contract  F30602-77-C-0163 ,  this  method  had  not  been  utilized  and  evalu¬ 
ated  for  any  conducting  geometries.  To  build  up  confidence  in  the  FD-TD  method 
for  future  applications,  it  was  desired  to  evaluate  the  usage  of  the  method 
for  certain  simple,  generic  metal  structures. 

During  the  research  program,  two  specific  metal  structures  were  used  in 
this  evaluation:  the  first,  an  aluminum  cylinder  with  one  open  end;  the 
second,  the  nose  cone  section  of  a  missile.  Each  structure  was  modeled  using 
the  FD-TD  method  to  compute  the  internal  EM  fields  generated  by  an  incident 
plane  wave  propagating  along  the  structure  axis.  The  results  of  the  cylinder 
model  were  then  compared  to  available  theoretical  and  experimental  data. 
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Final  evaluation  of  the  nose  cone  results  will  be  possible  when  reliable 
experimental  data  for  this  geometry  is  obtained  in  a  future  research  program. 

The  following  describes  each  coupling-analysis  problem  considered. 

2.3.2  Description  of  Coupling-Analysis  Problems  Considered 

Task  1:  Prediction  of  the  Coupling 
_ Into  an  Open-Ended  Cylinder 

The  FD-TD  technique  was  employed  to  solve  the  following  electromagnetic 
coupling  problem: 

Interacting  structure  -  Circular  (19.0  cm  diameter),  68.5  cm  long, 
open-ended  aluminum  cylinder,  as  shown  in  Figure  8  of  Final  Report 
RADC-TR-78-142  [Reference  8]; 

Incident  wave  -  300  MHz  plane  wave  propagating  down  the  cylinder  axis 
toward  its  open  end; 

Desired  fields  -  Each  component  of  total  F  and  total  IT  in  the  axial 
cross-section  plane  of  the  cylinder  down  to  40  cm  from  the  open  end. 
First,  with  the  cross-section  plane  parallel  to  the  incident  F,  and 
again  with  the  plane  parallel  to  the  incident  H; 

Resolution  -  0.5  cm  uniformly  throughout  the  mapping  planes; 

Plotted  values  -  In  decibels  relative  to  an  incident  F  of  1  volt/meter 
and  an  incident  H  of  1/377  ampere/meter. 

To  solve  this  coupling  problem,  an  existing  FD-TD  computer  code  was 
suitably  modified  to  model  the  cylinder.  Details  of  the  model  are  given  in 
Section  4.3.1  of  Reference  8. 

Task  2:  Prediction  of  the  Coupling 
_ Into  a  Missile  Nose  Cone 

The  FD-TD  technique  was  employed  to  solve  the  following  electromagnetic 
coupling  problem: 

Interacting  structure  -  Aluminum  nose  cone  (shown  in  Figure  9  of  Refer¬ 
ence  8)  with  two  apertures:  a  circular  one  in  the  nose,  and  a  sleeve 
fitting  located  23-1/3  cm  aft.  Missile  body  geometry  beyond  sleeve 
fitting  assumed  to  continue  to  infinity  with  constant  cross-section 
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shape.  Aperture  cases  investigated: 

Trial  1  -  Sleeve  fitting  open,  nose  aperture  closed, 

Trial  2  -  Both  apertures  open; 

Incident  wave  -  300  MHz  plane  wave  propagating  down  the  axis  of  the 
structure  toward  its  nose  aperture; 

Desired  fields  -  Each  component  of  total  E  and  total  FT  in  the  axial  cross- 
section  plane  of  the  nose  cone.  First,  with  the  cross-section  plane 
parallel  to  the  incident  E,  and  again  with  the  plane  parallel  to  the 
incident  H; 

Resolution  -  1/3  cm  uniformly  throughout  the  mapping  planes; 

Plotted  values  -  In  decibels  relative  to  an  incident  T  of  1  volt/meter 
and  an  incident  H  of  1/377  ampere/meter. 

To  solve  this  coupling  problem,  the  computer  code  of  Task  1  was  utilized 
with  a  new  set  of  data  cards  specifying  the  nose  cone  geometry.  Details  of 
the  model  are  given  in  Section  4.4.1  of  Reference  8. 

2.3.3  Results  and  Discussion 

Fields  computed  using  the  FD-TD  technique  closely  approximated  prior 
experimental  and  theoretical  data  (where  available)  for  the  cylinder  of 
Task  1.  Future  work  will  be  required  to  develop  comparative  experimental  data 
for  the  more  complex  nose  cone  geometry  of  Task  2.  Details  of  the  results  are 
given  in  Sections  4.3.3,  4.3.4,  and  4.4.2  of  Reference  8. 

This  research  program  demonstrated  that  the  FD-TD  method  can  be  success¬ 
fully  applied  to  axial -incidence,  electromagnetic  coupling  problems  involving 
highly  conducting  structures  with  hole  and  sleeve-type  apertures.  Accuracy 
of  the  FD-TD  results  was  very  good  relative  to  the  uncertainties  of  available 
experimental  and  numerical-theory  approaches.  Convergence  of  the  EM  fields 
to  the  sinusoidal  steady  state  occurred  within  about  2  cycles  of  the  incident 
wave.  This  resulted  in  program  running  times  of  3.5  minutes  or  less  on  the 
CDC  STAR-100  for  10^  -  cell  lattices. 

The  following  questions  concerning  the  FD-TD  method  remained  unanswered 
at  the  conclusion  of  RADC  Contract  F30602-77-C-0163: 
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1.  Can  the  FD-TD  method  accurately  model  EM  coupling  into  a  conducting 
structure  for  arbitrary  angles  of  incidence  and  arbitrary  wave 
polarization? 

2.  Can  the  FD-TD  method  accurately  model  a  complex  structure  with 
both  irregular-shaped  apertures  and  dielectric  or  permeable 
materials  within  the  interior? 

3.  Can  the  FD-TD  method  predict  the  voltage,  current,  or  power 
that  the  coupled  field  may  induce  on  wires  or  cable  bundles 
inside  of  a  metal  structure? 

4.  What  are  the  size  limits  of  an  object  that  can  be  modeled  using 
the  FD-TD  method? 

5.  What  is  the  capability  of  the  FD-TD  method  to  interface  with 
other  analysis  techniques,  such  as  "Method  of  Moments"? 

These  questions  should  be  addressed  systematically  to  permit  better  defini¬ 
tion  of  the  capabilities  and  limits  of  the  FD-TD  technique.  Favorable  con¬ 
clusions  on  most  or  all  of  these  questions  would  probably  lead  to  eventual 
application  of  the  FD-TD  method  to  complex  EM  coupling  problems  of  interest 
to  RADC  and  the  electromagnetic  compatibility  community. 

The  five  questions  above  probe  into  the  nature  of  the  FD-TD  modeling 
method  in  a  fundamental  manner.  For  example.  Question  1  involves  possibly 
modeling  the  longitudinal,  resonant  currents  set  up  in  an  elongated  metal 
structure  by  an  incident  TM  wave.  In  the  time  domain,  these  currents  may 
only  slowly  converge  to  the  sinusoidal  steady  state.  Answering  Question  1 
will  tell  us  how  much  time  is  needed  for  convergence,  and  what  measures,  if 
any,  can  be  used  to  speed  convergence  without  sacrificing  accuracy. 

Question  2  involves  the  interaction  of  conducting  and  dielectric  materi¬ 
als  subject  to  an  incident  wave.  In  the  time  domain,  the  presence  of  dielec¬ 
tric,  materials  within  a  conducting  enclosure  can  modify  the  convergence  of 
the  fields  within  the  enclosure  to  the  sinusoidal  steady  state.  Further, 
field  patterns  within  the  enclosure  may  be  altered  significantly  by  the  pres¬ 
ence  of  the  dielectric  media.  Answering  Question  2  will  tell  us  the  maximum 
level  of  complexity  that  the  FD-TD  method  can  attain  without  sacrificing  ac¬ 
curacy,  and  without  having  lengthy  running  times. 


26 


Question  3  involves  the  ability  of  the  FD-TD  method  to  model  conducting 
filaments  (wires)  instead  of  sheets.  Further,  this  question  will  probe  the 
modeling  of  wire  terminations  such  as  resistors,  and  the  ability  to  translate 
computed  field  quantities  such  as  E  and  H  to  circuit  quantities  such  as  V,  I, 
and  P. 

Question  4  involves  the  interaction  between  the  type  of  computer  being 
used  and  the  maximum  size  object  that  can  be  modeled  using  the  FD-TD  method. 
Here,  considerations  such  as  computer  cycle  time,  computer  magnetic-core 
memory  size,  computer  architecture  (pipeline  structure,  if  used;  word  size; 
memory  paging  procedure,  etc.),  and  overall  modeling  accuracy  desired  will 
be  considered.  This  study  will  tell  us  if  the  FD-TD  method  can  be  success¬ 
fully  and  efficiently  adapted  to  available  general  purpose  computers,  in 
addition  to  large-scale  processors  such  as  the  CDC  STAR  100. 

Lastly,  Question  5  involves  the  feasibility  of  developing  a  hybrid 
approach  where  some  appropriate  technique  (such  as  the  Method  of  Moments) 
is  used  to  develop  data  for  the  skin  currents  or  aperture  fields  of  a  large 
structure,  and  the  FD-TD  method  uses  this  data  to  compute  the  fields  pene¬ 
trating  small  apertures  of  the  structure  into  complex  interior  cavities.  In 
this  way,  each  technique  would  be  applied  in  the  structure-size/complexity 
range  that  it  is  best  suited  for,  allowing  an  overall  solution  that  is  ac¬ 
curate  for  large,  simple  structures  with  small,  complex  apertures  or  other 
fine  details.  Clearly,  adaptation  of  the  amplitude  and  phase  data  obtained 
from  a  frequency  domain  method  such  as  MOM  will  require  transformation  to 
the  time  domain.  Further,  usage  of  the  transformed  data  will  require  modif¬ 
ication  of  the  plane  wave  source  condition  normally  employed  for  the  FD-TD 
method. 

The  need  to  deal  with  the  five  questions  raised  above  forms  the  basis 
of  the  research  performed  under  the  present  contract,  F30602-79-C-0039.  The 
following  five  sections  will  discuss  the  results  of  the  present  research. 
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3.0  PROGRAM  TASK  1 :  PREDICTION  OF  THE  COUPLING 
_ INTO  AN  OPEN-ENDED  CYLINDER 

3.1  Broadside  Incidence  and  Transverse  Electric  (TE) 

Wave  Polarization _ 

The  FD-TD  technique  was  employed  to  solve  the  following  electromagnetic 
coupling  problem: 

Interacting  structure  -  Circular  (19.0  cm  diameter),  68.5  cm  long, 
open-ended  aluminum  cylinder  as  shown  in  Figures  3  and  4; 

Incident  wave  -  300  MHz  plane  wave  propagating  perpendicular  to  the  cyl¬ 
inder  axis  and  polarized  so  that  the  wave  electric  field  is  perpendicular 
to  the  plane  formed  by  the  cylinder  axis  and  the  incident  wavevector; 

Desired  fields  -  Each  component  of  total  E  and  total  H  in  the  axial 
cross-section  plane  of  the  cylinder  down  to  40  cm  from  the  open  end. 
First,  with  the  cross-section  plane  parallel  to  the  incident  E\  and 
again  with  the  plane  parallel  to  the  incident  wavevector; 

Resolution  -  0.5  cm  uniformly  throughout  the  mapping  planes; 

Plotted  values  -  In  decibels  relative  to  an  incident  E  of  1  volt/meter 
and  an  incident  H  of  1/377  ampere/meter. 

To  solve  this  coupling  problem,  an  existing  FD-TD  computer  code  was 
suitably  modified  to  model  the  cylinder  of  Figure  3  and  4.  The  following 
steps  were  taken: 

a.  A  160  x  63  x  24  cell  lattice  (1.45-106  unknown  field  components) 
was  programmed,  with  even  symmetry  of  the  incident  fields  and 
cylinder  assumed  about  lattice  plane  z  =  246. 

b.  The  cylinder  geometry  of  Figures  3  and  4  was  mapped  into  the 
finite-difference  lattice  for  a  unit  cell  diameter  of  6  =  0.5  cm. 

c.  The  FD-TD  program  was  run  for  this  geometry  for  800  time  steps 
(equivalent  to  2.0  cycles  of  the  incident  wave)  assuming  a  slight 
amount  of  isotropic  loss  (o.^  =  0.01  mho/m)  for  the  air  within 

the  cylinder  to  increase  the  rate  of  convergence  to  the  steady  state. 
2.5*1 0^  words  of  memory  and  9.4  minutes  of  central  processor  time 
were  required  on  the  Control  Data  STAR-100. 
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Circular  cylinder 
cross-section 


Fig.  4  CROSS-SECTION  OF  MODEL  OF  OPEN-ENDED 
ALUMINUM  CYLINDER,  TE/BROADSIDE  CASE 
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;  d.  The  computed  results  for  E  and  H  were  reduced  to  contour  maps  along 

the  horizontal  symmetry  plane  and  the  vertical  observation  plane. 

Figure  5  plots  the  computed  total  electric  field,  |Ez/Ez  I  in  decibels 
\  along  the  cylinder  axis  for  the  cases  n  =  200,  n  =  400,  n  =  6^  and  n  =  800 

time  steps.  Each  curve  gives  the  computed  field  envelope  during  the  200  time- 

Sstep  period  (0.5  cycle  of  the  incident  wave)  before  the  specified  value  of  n. 

In  Figure  5,  it  should  be  noted  that,  at  200  time  steps,  the  incident  wave 
penetrated  only  about  40  cm  into  the  cylinder,  causing  the  sharp  downward 
break  of  the  n  =  200  curve.  For  all  other  curves,  the  incident  wave  penetrat¬ 
ed  fully  to  the  cylinder  backplane. 

Figure  5  indicates  the  rate  of  field  convergence  to  the  sinusoidal  steady 

state  under  the  condition  of  finite  a.  ..  Once  the  wave  fields  were  estab- 

i  nt 

lished  throughout  the  cylinder  (curves  n  =  400,  600,  and  800),  the  principal 
effect  of  an  added  200  program  time  steps  was  the  lengthening  of  the  approx¬ 
imately  linear  decibel  slope  extending  from  the  aperture,  and  a  consequent 
deepening  of  the  field  null  within  the  cylinder.  Noting  the  break  points  of 
the  curves  from  the  linear  decibel  slope,  each  200  time-step  increase  is  seen 
to  have  lengthened  the  slope  by  about  10  cm  and  decreased  the  residual  inter- 
j  nal  fields  by  about  10  dB. 

|  Figure  6  compares  the  FD-TD  computed  |Hx/H  |  along  the  cylinder  axis 

$  with  computed  results  by  D.  Wilton  and  A.  Glisson  using  a  frequency-domain, 

body-of-revolution  computer  program  [9].  The  FD-TD  results,  shown  as  a  solid 
curve,  are  after  800  time  steps.  The  frequency-domain-program  results  are 
graphed  as  a  dashed  curve.  Excellent  agreement  (within  0.5  dB)  of  the  results 
of  the  two  approaches  is  seen  over  the  first  15  cm  of  field  penetration  into 
the  cylinder.  Over  this  span,  the  total  Hx  field  decay  is  of  the  order  of 

50-55  dB.  The  rate  of  decay  of  H  computed  by  either  technique  equals  3.29 

/\ 

dB/cm,  which  compares  favorably  with  the  3.46  dB/cm  rate  predicted  by  simple 
waveguide  mode  theory  for  this  beyond-cutoff  case. 

Using  results  of  the  FD-TD  program  (800  time  steps,  a^nt  =  0.01  mho/m), 
Figures  7-9  graph  contour  maps  of  the  three  computed  field  components  at  the 
cylinder's  horizontal  symmetry  plane.  Figures  10  and  11  graph  contour  maps 
of  the  six  computed  field  components  at  the  cylinder's  vertical  observation 
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CONVERGENCE  OF  THE  COMPUTED  ELECTRIC  FIELD 
ALONG  THE  CYLINDER  AXIS,  TE/BROADSIDE  CASE 
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COMPARISON  OF  RESULTS  FOR  THE  AXIAL  MAGNETIC  FIELD 
WITHIN  THE  CYLINDER,  TE/ BROADSIDE  CASE 
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Fig. 10.  COMPUTED  ELECTRIC  FIELD  CONTOURS  IN  THE 

CYLINDER  VERT.  OBS .  PLANE,  TE/ BROADSIDE  CASE 
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plane.  For  optimum  clarity,  contours  are  plotted  at  either  6  dB  or  10  dB 
intervals  in  these  figures.  Contour  points  were  generated  by  using  a  linear 
interpolation  method  to  determine  their  locations  between  adjacent  FD-TD 
field  envelope  points.  Although  the  lattice  cell  size,  6,  equals  0.5  cm, 
this  interpolation  method  allows  the  generation  of  smooth  curves  in  most 
cases  without  a  0.5  cm-period  stair-case  effect. 


3.2  Broadside  Incidence  and  Transverse  Magnetic  (TM) 

Wave  Polarization _ 

The  FD-TD  technique  was  employed  to  model  the  same  broadside-incidence 
cylinder  problem  as  that  of  Section  3.1  and  Figures  3  and  4,  except  that  the 
incident  wave  polarization  was  changed  to  align  the  wave  magnetic  field  per¬ 
pendicular  to  the  plane  formed  by  the  cylinder  axis  and  the  incident  wave- 
vector.  To  solve  this  coupling  problem,  the  FD-TD  computer  code  of  Section 
3.1  was  modified  to  accept  the  polarization  change.  Further,  the  code  was 
improved  to  reduce  central  processing  (CP)  time  by  about  48%.  This  was  ac¬ 
complished  by  storing  in  memory  an  array  of  constant  coefficients  which  pre¬ 
viously  had  been  re-computed  each  time  step.  For  example,  the  TM  case  re¬ 
quired  only  4.9  minutes  of  STAR-100  CP  time  for  800  time  steps  of  1.45-10® 
unknown  field  components,  as  compared  to  9.4  minutes  of  CP  time  for  the  cor¬ 
responding  TE  case  of  Section  3.1.  All  FD-TD  programs  discussed  subsequently 
in  this  report  incorporate  the  improved  code. 

Figure  12  plots  the  FD-TD  computed  total  radial  magnetic  field, 

|H  /H  |  (in  decibels),  along  the  cylinder  axis  for  the  cases  n  =  200, 
z  zinc 

n  =  400,  n  =  600,  and  n  =  800  time  steps.  As  such.  Figure  12  indicates  the 
rate  of  field  convergence  to  the  sinusoidal  steady  state.  Comparing  Figure 
12  to  Figure  5,  it  is  seen  that  the  internal  field  for  the  TM  case  converges 
more  rapidly  to  the  steady  state  than  the  field  for  the  TE  case.  This  is  a 
somewhat  surprising  result  since  the  TM  polarization  is  expected  to  induce 
damped,  oscillating,  longitudinal  current  flow  along  the  cylinder  which  might 
be  expected  to  delay  convergence  to  the  steady  state.  Evidently,  the  partic¬ 
ular  cylinder  electrical  length  and  effective  radiation  resistance  at  300  MHz 
combined  to  cause  rapid  damping  of  these  oscillating  currents,  so  ihat  the 


39 


n  =  400  curve 


n  * 


i  - - - ■+ - 1  ■>  ■  —  < - -H 

L0  20  30  40  50  60 

Distance  from  aperture  (cm) 


CONVERGENCE  OF  THE  COMPUTED  MAGNETIC  FIELD 
ALONG  THE  CYLINDER  AXIS,  TM/BROADSIDE  CASE 


"W 


internal  fields  could  rapidly  converge. 

Figure  13  compares  the  FD-TD  computed  total  axial  electric  field, 

! E  / E  |,  along  the  cylinder  axis  with  computed  results  by  D.  Wilton  and 
x  xinc 

A.  Glisson  using  a  frequency-domain,  moment-method,  body-of-revolution  com¬ 
puter  program  [9].  The  FD-TD  results,  shown  as  a  solid  curve,  are  after  800 
time  steps.  The  frequency-domain-program  results  are  graphed  as  a  dashed 
curve.  Except  near  the  aperture,  good  agreement  (within  2.5  dB)  of  the 
results  of  the  two  approaches  is  seen  over  the  first  25  cm  of  field  penetra¬ 
tion  into  the  cylinder.  Over  this  span,  the  total  Ex  field  decay  is  of  the 
order  of  50-55  dB.  The  rate  of  decay  of  computed  by  either  technique 
equals  2.13  dB/cm,  which  is  exactly  the  rate  predicted  by  waveguide  mode 
theory  for  this  beyond-cutoff  case. 

Using  results  of  the  FD-TD  program  (800  time  steps,  o-nt  =  0.01  mho/m), 
Figures  14-16  graph  contour  maps  of  the  three  computed  field  components  at  the 
cylinder's  horizontal  symmetry  plane.  Figures  17  and  18  graph  contour  maps 
of  the  six  computed  field  components  at  the  cylinder's  vertical  observation 
plane.  For  optimum  clarity,  contours  are  plotted  at  either  6  dB  or  10  dB  in¬ 
tervals  in  these  figures. 

3.3  Oblique  Incidence  and  Transverse  Magnetic  (TM)  Wave  Polarization 

As  discussed  in  Section  7.0  of  this  report,  it  appears  possible  to  ac¬ 
curately  model  aperture-excited  structures  via  a  hybrid  method-of-moments/ 
FD-TD  approach.  This  hybrid  approach  has  the  advantage  of  allowing  arbitrary 
field  excitation  without  requiring  re-programming  of  the  basic  geometry  of 
a  particular  interacting  structure.  Further,  it  allows  small,  highly-complex 
interior  zones  of  large,  simple  structures  to  be  analyzed. 

As  one  of  the  test  cases  for  the  hybrid  method,  it  was  decided  to  model 
the  field  penetration  within  the  open-ended  cylinder  of  Section  3.1  and  3.2, 
but  assuming  that  the  incident  wave  propagates  at  an  angle  of  45°  to  the  cyl¬ 
inder  axis  toward  the  open  end.  The  wave  was  assumed  to  be  polarized  trans¬ 
verse  magnetic  (TM),  i.e.,  the  wave  magnetic  field  was  perpendicular  to  the 
plane  formed  by  the  cylinder  axis  and  the  incident  wavevector.  Results  for 
this  case  are  documented  in  Section  7.4.  These  results  show  a  very  good 
agreement  between  the  hybrid  method  and  conventional  moment-method  results, 
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(decibels) 


Fig.  13.  COMPARISON  OF  RESULTS  FOR  THE  AXIAL  ELECTRIC  FIELD 
WITHIN  THE  CYLINDER,  TM/ BROADSIDE  CASE 
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Fig.  16.  COMPUTED  E  CONTOURS  IN  CYLINDER  HOR. 
SYMMETRY  PLANE,  TM/ BROADS IDE  CASE 
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Fig.  18. 


COMPUTED  MAGNETIC  FIELD  CONTOURS  IN  THE 
CYLINDER  VERT.  OBS .  PLANE,  TM/BROADSIDE  CASE 


p 


indicating  the  feasibility  of  modeling  structures  for  arbitrary  angles  of 
wave  incidence. 

3.4  Conclusions 

Completion  of  Program  Task  1  demonstrates  the  capability  of  the  FD-TD 
method  to  accurately  model  the  electromagnetic  fields  coupled  into  a  structure 
which  is  illuminated  with  an  electromagnetic  plane  wave  having  an  arbitrary 
polarization.  Elements  of  Task  5  will  demonstrate  this  modeling  capability 
for  an  arbitrary  wave  angle  of  incidence,  as  well. 
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4.0  PROGRAM  TASK  2:  PREDICTION  OF  THE  COUPLING  INTO  A  MISSILE 

GUIDANCE  SECTION 


4 . I  Model  1  -  Includes  the  Interior  Dielectric  Material 

The  FD-TD  technique  was  employed  to  solve  the  following  electromagnetic 
coupling  problem: 

Interacting  structure  -  Empty  shell  of  missile  guidance  section,  shown 
in  Figures  19-21,  with  two  apertures:  a  circular  one  in  the  nose,  and  a 
sleeve  fitting  located  23-1/3  cm  aft.  Missile  body  geometry  beyond  sleeve 
fitting  assumed  to  continue  to  infinity  with  constant  cross-section  shape. 
Dielectric  materials  are  located  as  shown  in  Figures  19  and  20.  Assumed 
electrical  parameters  for  the  media  comprising  the  model: 


Med i urn 

Rel.  Permittivity,  £  r 

Conductivity,  a  (mhos/m) 

Aluminum 

1.0 

3.7  -107 

Fiberglas 

5.5 

2.4  -10‘3 

Phenol ic 

4.5 

1 

o 

o 

CO 

Magnesium  fluoride 

5.3 

0.0 

For  the  fiberglas,  phenolic,  and  magnesium  fluoride,  the  effective  conductiv¬ 
ity,  a,  is  given  by 

a  =  use"  =  (  2tt  •  3  •  10^)  •  (e^  •  cq  •  tan  5) 

=  0.0167  E.J,  tanS  mhos/m 
where  tan  6  is  the  loss  tangent. 

Incident  wave  -  300  MHz  plane  wave  propagating  down  the  axis  of  the  structure 
toward  its  nose  aperture; 

Desired  fields  -  Each  component  of  total  E  and  total  H  in  two  axial  cross- 
section  planes  of  the  structure.  First,  with  the  cross-section  plane  parallel 
to  the  incident  E,  and  again  with  the  plane  parallel  to  the  incident  H; 

Resolution  -  1/3  cm  uniformly  throughout  the  mapping  planes; 

Plotted  values  -  In  decibels  relative  to  an  incident  E  of  1  volt/meter  and  an 
incident  H  of  1/377  ampere/meter. 

To  solve  this  coupling  problem,  the  following  steps  were  taken: 
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aluminum  backplane 


Fig. 19.  GEOMETRY  OF  NOSE  CONE  MODEL  AT  VERTICAL  SYMMETRY  PLANE, 
SHOWING  COMPONENT  MATERIALS. 
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Fig. 20  .  GEOMETRY  OF  NOSE  CONE  MODEL  AT  HORIZONTAL  SYMMETRY  PLANE, 
SHOWING  COMPONENT  MATERIALS. 
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Fig.  21.  KEY  CROSS-SECTIONS  OF  NOSE-CONE  MODEL 


a.  The  FD-TD  computer  program  modeling  the  metal-only  guidance 
section  shell  of  RADC  Contract  F30602-77-C-0163  [Reference  8] 
was  re-written  to  incorporate  the  algorithm  speed-up  provisions 
discussed  in  Section  3.2. 

b.  A  new  data  card  deck  (approximately  700  cards)  was  constructed 
to  specify  the  metal  and  dielectric  structure  of  Figures  19-21. 

c.  The  FD-TD  program  was  run  for  1800  time  steps  (equivalent  to  3.0 
cycles  of  the  incident  wave),  assuming  a  slight  amount  of  isotropic 
loss  (ck  =  0.025  mho/m)  for  the  air  within  the  cylinder  to  speed 
the  rate  of  convergence.  7.8  x  10^  words  of  memory  and  3.7  minutes 
of  STAR-100  central  processor  time  were  required. 

d.  The  computed  results  for  E  and  H  were  reduced  to  contour  maps  along 
the  symmetry  planes. 

Figure  22  plots  the  computed  contours  for  the  3  field  components  in  the 
vertical  symmetry  plane,  while  Figure  23  plots  the  contours  for  the  3  field 
components  in  the  horizontal  symmetry  plane.  These  figures  can  be  directly 
compared  with  Figures  23  and  24,  respectively,  of  Reference  8,  which  graph 
the  FD-TD  computed  field  contours  for  the  missile  guidance  section  shell  with 
out  dielectrics.  An  intriguing  observation  is  that  the  presence  of  the  fiber 
glas  in  the  circumferential  sleeve  fitting  serves  to  increase  the  interior 
field  level  toward  the  rear  of  the  guidance  section  shell  by  an  average  of 
about  15  dB.  The  general  shape  of  the  contours  is  unchanged,  however. 

This  increase  at  first  seems  to  be  anomalous,  since  one  would  ordinarily 
expect  that  loading  an  aperture  with  a  high-permittivity  dielectric  such  as 
fiberglas  would  tend  to  reduce  the  transmitted  fields.  Flowever,  two  reason¬ 
able  physical  mechanisms  could  explain  the  apparent  anomaly: 

1.  The  fiberglas-loaded  sleeve  fitting  may  form  a  nearly  resonant 
aperture  at  300  MHz.  This  is  because  the  wavelength  in  fiberglas 
at  this  frequency  is  42.6  cm  (assuming  a  fiberglas  relative  permit¬ 
tivity  of  5.5),  very  nearly  equal  to  the  40.2-cm  physical  circum¬ 
ference  of  the  sleeve  fitting.  Thus,  any  radial  excitation  of  the 
sleeve-fitting  gap  would  tend  to  be  reinforced  in  phase  along  the 
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Fig.  22.  COMPUTED  FIELD  CONTOURS  IN  NOSE  CONE  VERTICAL  SYMMETRY  PLANE 
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23.  COMPUTED  FIELD  CONTOURS  IN  NOSE  CONE  HORIZONTAL 
SYMMETRY  PLANE. 


circumference  of  the  sleeve  fitting. 

2.  The  fields  excited  radially  in  the  fiberglas-loaded  sleeve 
fitting  may  be  channeled  into  the  interior  of  the  guidance 
section  via  the  fiberglas  itself,  which  might  form  a  surface 
wave-guiding  action  in  conjunction  with  the  metal  surrounding 
it. 

The  possible  existence  of  the  resonant-aperture  phenomenon  could  be  tested 
numerically  by  re-running  the  FD-TD  model  for  several  frequencies  above  and 
below  300  MHz.  If  the  resonance  exists,  field  penetration  should  be  maximum 
at  one  frequency  and  fall  off  to  either  side. 

4.2  Model  2  -  Includes  the  Key  Interior  Metal  Components  as  Well 
as  the  Dielectric  Material _ 

The  FD-TD  technique  was  employed  to  solve  the  following  electromagnetic 
coupling  problem: 

Interacting  structure  -  Identical  to  that  of  Model  1  (Figures  19-21),  with 
the  addition  of  the  following  metal  and  dielectric  interior  components: 

1)  Head  coil  assembly  (assumed  solid  metal); 

2)  Cooled  detector  unit  -  CDU  (assumed  solid  metal); 

3)  Phenolic  ring  around  the  CDU; 

4)  Preamp  can  (metal ) ; 

5)  Wire  connecting  the  CDU  to  the  preamp  can; 

6)  Wire  connecting  the  preamp  can  to  the  metal  backplane;  and 

7)  Longitudinal  metal  support  rods. 

These  additional  components  are  located  as  shown  in  Figures  24-26. 

Incident  wave  -  Axially  directed  300  MHz,  as  in  Model  1; 

Desired  fields  -  Total  E”  and  total  H  in  two  axial  cross-section  planes, 
as  in  Model  1 . 
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Fig. 25.  GEOMETRY  OF  NOSE  CONE  MODEL  AT  HORIZONTAL  OBSERVATION  PLANE, 
SHOWING  COMPONENT  MATERIALS. 


Fig. 26.  KEY  CROSS-SECTLONS  OF  NOSE  CONE  MODEL.  HEAD-COIL  ASSEMBLY  AND 
COOLED  DETECTOR  UNIT  NOT  SHOWN. 


Resolution  -  1/3  cm,  as  in  Model  1. 


Plotted  values  -  In  decibels  relative  to  incident  E  and  incident  H,  as 
in  Model  1 . 

To  solve  this  coupling  problem,  the  following  steps  were  taken: 

a.  The  FD-TD  computer  program  for  Model  1  of  the  missile  guidance 
section  was  re-written  to  eliminate  the  horizontal  symmetry  plane, 
since  the  interior  metal  components  do  not  have  this  symmetry. 

b.  Drawings  of  the  missile  guidance  section  were  studied  to  determine 
the  interior  geometry.  A  new  data  card  deck  (approximately  700 
cards)  was  constructed  to  specify  the  geometry. 

c.  The  FD-TD  program  was  run  for  1800  time  steps  (equivalent  to  3.0 
cycles  of  the  incident  wave),  assuming  a  slight  amount  of  isotropic 
loss  (a.  t  =  0.025  mho/m)  for  the  air  within  the  cylinder  to  speed 
the  rate  of  convergence.  1.6 ‘10®  words  of  memory  and  7.0  minutes 
of  STAR-100  central  processor  time  were  required. 

d.  The  computed  results  for  F  and  H  were  reduced  to  contour  maps  along 
the  symmetry  planes. 

Figures  27-29  plot  the  computed  contours  for  the  3  field  components  in 
the  vertical  symmetry  plane,  while  Figures  30  and  31  plot  the  contours  for 
the  6  field  components  in  the  horizontal  observation  plane.  These  figures 
can  be  directly  compared  with  Figures  22  and  23  for  Model  1.  It  is  seen  that 
the  metal  components  and  wires  serve  to  significantly  disturb  the  field  con¬ 
tours  from  the  empty-shell  case.  Further,  greatly  increased  fields  penetrate 
into  the  interior  of  the  guidance  section  near  the  metal  components  and  wires. 
An  important  observation  is  that  the  wires  connecting  the  cooled  detector 
unit,  preamp  can,  and  metal  backplane  are  paralleled  by  high-level  magnetic 
field  contours  (as  evidenced  in  Figures  28,  31a,  and  31b).  This  is  indica¬ 
tive  of  substantial,  uniform  current  flow  along  these  wires,  and  will  be  dis¬ 
cussed  in  more  detail  in  the  next  section. 
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Fig. 28.  COMPUTED  H  CONTOURS  IN  NOSE  CONE  VERTICAL  SYMMETRY  PLANE  (IN  DB) 
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Fig.  31.  MAGNETIC  FIELD  CONTOURS  IN  NOSE  CONE  HORIZONTAL  OBS .  PLANE  (IN  DB) 
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4.3  Conclusions 

Completion  of  Program  Task  2  demonstrates  the  capability  of  the  FD-TD 
method  to  map  fields  coupled  into  a  complex  structure  that  has  both  irregular 
shaped  apertures  and  interior  dielectric  and  metallic  materials.  Determina¬ 
tion  of  the  accuracy  bound  for  this  method,  as  applied  to  complex  structures, 
awaits  the  resu'ts  of  future  experimental  programs  since  other  numerical  ap¬ 
proaches  cannot  deal  directly  with  this  level  of  complexity.  No  comparative 
data  is  available  at  this  time. 
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5.0  PROGRAM  TASK  3:  PREDICTION  OF  THE  COUPLING  TO  WIRES 

OR  CABLE  BUNDLES 


A  study  of  the  ability  of  the  FD-TD  method  to  predict  electromagnetic 
coupling  to  a  wire  or  cable  bundle  was  made  using  Model  2  of  the  missile 
guidance  section  (detailed  in  Section  4.2).  As  stated  in  that  section  of 
the  report,  the  model  included  two  connecting  wires:  1)  Wire  1  -  between 
the  cooled  detector  unit  and  the  preamp  can;  and  2)  Wire  2  -  between  the 
preamp  can  and  the  metal  backplane.  These  wires  were  really  idealizations  of 
a  more  complex  situation  in  which  two  multi  conductor  wire  bundles  extended 
between  the  structures  mentioned.  Some  wires  of  the  bundles  terminated  in 
electronic  components;  some  were  small-diameter  cables  with  grounded  shields. 
The  goal  was  to  model  the  direct  metal -to-metal  current  path  through  the  mis¬ 
sile  guidance  section.  Therefore,  using  a  simple  single-wire  model  for  each 
bundle,  only  the  ground-path  current  and  equivalent  common-mode  current  of 
each  bundle  was  modeled  using  the  FD-TD  method. 

5. 1  Wire  Model 

Modeling  of  a  connecting  wire  using  the  FD-TD  method  can  be  done  best  in 
a  direct  manner  by  specifying  the  wire  model  as  a  filamentary  conducting  path 
of  realistic  location  relative  to  other  components  of  the  model.  In  this 
manner,  certain  transmission  line  properties  of  the  wire,  such  as  capacitance 
and  external  inductance,  that  arise  due  to  its  location  relative  to  other 
structures  are  automatical ly  accounted  for  in  the  time-domain  solution  of 
Maxwell's  equations  for  the  whole  model,  which  includes  the  interaction  of 
the  wire  and  all  adjacent  structures.  The  transmission-line  properties  of 
the  wire  that  arise  due  to  internal  impedance  are  accounted  for  by  adjusting 
the  model  wire  conductivity  so  that  these  properties  match  those  of  the 
actual  conductor.  This  procedure  is  now  discussed  for  the  wire  model  of 
Section  4.2. 

For  the  single-wire  equivalent  of  each  wire  bundle,  it  was  decided  to 
select  a  round  copper  wire  having  an  American  (B  &  S)  gauge  of  15.  The 
internal  impedance  of  this  wire  at  high  frequencies  [10]  was 
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09a) 


I 


Zi. 


=  Rs  (1  +  j) 


wire 


2nr 


=  (2.61 *10-Z  /T)  •  (1  +  j) 
(2tt)  •  (0.000725) 


=  0.992  •  (1  +  j)  ohms/meter 
where  Rg  =  wire's  copper  surface  resistivity  (in  ohms); 
f  =  300  MHz; 

rQ  =  wire  radius  =  0.000725  meters. 

For  this  choice,  the  internal  resistance  at  300  MHz  was 

R.j  =1.0  ohm/meter 

wire 

and  the  internal  reactance  was 

X.  =  +1.0  ohm/meter, 

wire 

For  the  FD-TD  model  of  the  above  wire,  a  conducting  filament  having  a 
square  cross  section  defined  by  a  single  lattice  unit  cell  (1/3  cm)  was  used 
By  choosing  the  conductivity  of  the  model  wire  medium  to  equal  a  =  0.66  •  10 
mhos/meter,  the  model  wire  could  be  defined  to  have  the  same  internal  imped¬ 
ance  as  the  above  15  gauge  copper  wire: 

Z,  =  R  •  (1  +  J)  (19b) 

model  model _ 

(model  wire  circumference) 


Trfuo 


model  wire 


(1  +  0) 


(2#) 


=  1.0  *  (1  +  j)  ohms/meter 

where  =  vacuum  permeability  =  47r-10’Z  henrys/meter 
f  =  300  MHz 

°model  wire  '  °-66'1C|7  ",hos/l"eter- 
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5.2  Results 


In  Figure  28  of  Section  4.2,  note  how  the  computed  0  dB  contour  of  H 
nearly  parallels  the  model  ground  wire  connecting  the  cooled  detector  unit  to 
the  pre-amp  can.  This  behavior,  combined  with  computed  data  showing  that  the 
longitudinal  electric  field  at  the  wire  is  small,  is  indicative  of  a  trans¬ 
verse  electromagnetic  (TEM)  wave  mode  propagating  along  the  ground  wire.  In 
effect,  current  flow  through  the  wire  generates  locally  a  magnetic  field  loop¬ 
ing  around  the  wire,  which  when  "cut"  by  the  vertical  symmetry  plane,  shows  up 
as  parallel  field  contours  spaced  equally  on  each  side  of  the  wire.  For  this 
case.  Ampere's  generalized  law  reduces  to  the  situation  illustrated  in  Figure 
32: 


(20a) 


(20b) 


In  this  manner,  Iwire  can  be  estimated  by  knowledge  of  the  magnitude  and  loca¬ 
tion  of  the  parallel,  adjacent,  magnetic  field  contours.  This  applies  equally 
well  to  the  ground  wire  connecting  the  pre-amp  can  to  the  backplane,  which  is 
seen  to  have  parallel  +  10  dB  Hx  contours  on  either  side. 

Table  1  lists  the  predicted  values  of  1^  (current  in  the  wire  from  the 

cooled  detector  unit  to  the  pre-amp  can)  and  ^  (current  in  the  wire  from  the 

pre-amp  can  to  the  backplane).  For  the  case  of  1^,  it  is  assumed  that 

H  =  H.  „  and  r  =  1.5  6  =  0.005  m.  For  the  case  of  I0  it  is  assumed  that 
me  o  2, 

H  =  3.16  H.  and  r  =1.55=  0.005  m. 
me  o 

Table  1 

PREDICTED  GROUND  WIRE  CURRENTS 


Einc  <v/0,) 

Hi„c  <A/"> 

PinC  (uW/ cm3) 

(uA) 

I2  (pA) 

0.1 

2.65-10"4 

2.65-10'3 

8.3 

26.3 

1  .0 

2.65' 10-3 

0.265 

83 

263 

10.0 

2.65' 10~3 

26.5 

830 

2630 

Current-carrying  wire 


Fig.  32.  GEOMETRY  OF  AMPERE'S  LAW 


5.3  Conclusions 


Completion  of  Program  Task  3  demonstrates  the  ability  of  the  FD-TD  tech¬ 
nique  to  directly  predict  the  current  that  the  coupled  fields  may  induce  on  a 
wire  or  bundle  of  wires  which  is  located  inside  a  complicated  structure  (such 
as  a  missile  guidance  section).  Here,  the  observed  parallelism  of  the  com¬ 
puted  magnetic  field  contours  near  the  modeled  wire(s)  is  employed  in  combina¬ 
tion  with  Ampere's  law  to  translate  FD-TD-computed  magnetic  fields  to  wire 
current  magnitudes.  Confirmation  of  the  FD-TD  computed  wire  currents  is  re¬ 
quired  in  a  future  experimental  program. 
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6.0  PROGRAM  TASK  4:  DETERMINATION  OF  THE  SIZE  LIMITS  OF  AN  OBJECT 

_ THAT  CAN  BE  MODELED  BY  THE  FD-TD  METHOD 

6. 1  Generdl  Considerations 

Like  many  other  electromagnetic-analysis  computer  methods,  the  FD-TD 
approach  requires  significant  quantities  of  computer  memory  and  processing 
time  to  model  realistic  problems.  In  general,  the  computer  resources  needed 
by  the  FD-TD  method  increase  for  larger  structures  and  finer  spatial  resolu¬ 
tions,  since  more  memory  and  processing  time  is  required  to  account  for  high¬ 
er  levels  of  geometric  detail.  Clearly,  arbitrarily-large  objects  cannot  be 
modeled  using  this  method. 

The  purpose  of  this  section  is  to  allow  an  estimate  of  the  maximum  size 
of  an  object  that  can  be  modeled  using  the  FD-TD  method,  given  limitations 
of  the  available  computer  system.  It  will  be  shown  that  the  most  important 
limiting  factors  of  any  computer  system  in  implementing  the  FD-TD  method  are: 

1.  Quantity  of  available  fast-access  central  memory,  such  as 
magnetic  core  and  semiconductor  arrays; 

2.  Data  transmission  speed  of  on-line  mass  storage  memory, 
such  as  magnetic  disc;  and 

3.  Floating-point  operation  speed. 

6.2  Required  Computer  Memory 

At  each  time  step  of  the  FD-TD  algorithm,  the  following  twelve  data 
quantities  are  stored  in  computer  memory  for  each  unit  cell  of  the  space  lat¬ 
tice:  1)  Three  vector  components  of  IT;  2)  Three  vector  components  of  H; 
and  3)  Six  constant  coefficients.  Additional  computer  memory  is  needed  to 
implement  the  lattice  truncation  conditions  at  the  outermost  cells  of  the 
lattice.  Overall,  it  has  been  found  that  the  following  is  an  accurate  estim¬ 
ate  of  the  total  computer  memory  requirement,  M,  for  a  FD-TD  lattice  contain¬ 
ing  N  cells: 

M  =  13.6N  words  of  memory.  (21a) 

The  N  cells  are  arranged  in  a  three-dimensional  cubic  lattice  to  permit 
the  maximum  flexibility  in  specifying  an  arbitrary  object.  Denoting  the  max- 
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imum  span  of  the  object  in  three  perpendicular  directions  as  D-^  ,  D^, 
meters,  respectively;  the  lattice  cell  size  as  6  meters;  and  the  number  of 
perpendicular  symmetry  planes  of  the  object/electromaqnetic  excitation  as 
s  =  0,  1 ,  or  2 ,  M  is  computed  to  be 

M  =  13.6N 

-  13.6  x  (^1)  x  (°2l  x  (^3) 

5 '  6 

=  13.6  D-,  D2  D3 

-  words  of  memory.  (21b) 

2s  6 3 

Equation  21b  specifies  the  memory  requirement  for  plane  wave  excitations 
having  arbitrary  polarization  and  angle  of  incidence.  Further,  using  the 
hybrid  computation  approach  developed  in  Section  7  of  this  report.  Equation 
21b  holds  for  both  compact,  isolated  objects  as  well  as  sub-zones  or  cavities 
of  extended  objects  which  can  be  treated  via  the  aperture  electric  current 
equivalence  theorem.  In  the  former  case,  D-j ,  D^,  and  are  simply  the  max¬ 
imum  spans  of  the  isolated  object.  In  the  latter  case,  D-j ,  D^,  and  are 
the  maximum  spans  of  only  the  sub-zone  or  cavity,  rather  than  the  object  as 
a  whole.  Thus,  the  hybrid  approach  allows  excellent  resolution  of  distinct 
zones  or  cavities  of  objects  that  may  have  overall  much  larger  dimensions,  as 
will  be  shown  in  Section  7.0. 

6.3  Required  Computer  Processing  Time 

6.3.1  Central  Processing  Time 

At  each  time  step  of  the  FD-TD  algorithm,  the  following  arithmetic 
operations  are  performed  for  each  unit  cell  of  the  space  lattice  (assuming 
permeability,  u,  equal  to  uQ  everywhere  in  the  lattice):  1)  6  floating¬ 
point  multiplications;  and  2)  24  floating-point  additions  or  subtractions. 

Additional  operations  are  needed  to  implement  the  lattice  truncation  condi¬ 
tions  and  the  envelope  computations.  Overall,  however,  it  has  been  found  that 
the  central  processing  time  is  dominated  by  the  basic  time-stepping  opera¬ 
tions. 
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Scalar  Processing  Machines 

Denoting  the  computation  time  for  one  single-precision  floating-point 
multiplication  as  x  ,  and  the  computation  time  for  one  single-precision 
floating-point  addition  as  t  ,  the  following  is  an  estimate  of  the  central 

a 

processing  time,  CP^ ,  for  one  time  step  of  the  FD-TD  algorithm  over  N  cells: 

CP,  -  (6  t_  +  24  t_ )  x  N  seconds  .  (22a) 

1  m  a  time  step 


Assuming  that  R  cycles  of  the  excitation  frequency,  f,  must  be  time-stepped 
to  attain  the  sinusoidal  steady  state,  the  number  of  FD-TD  algorithm  time 
steps  needed  is 


nmax  -  (^)  _  (^r)  _  ^  time  steps 


(22b) 


6t 


(A) 

2C  ‘ 


where  6t  is  the  time-step  size  of  the  FD-TD  algorithm,  6  is  the  lattice-cell 
size,  and  c  is  the  speed  of  light  in  vacuum.  Therefore,  from  Equations  22a 
and  22b,  the  total  central  processing  time,  CP,  is  estimated  as 

CP  =  CP,  x  n 

1  max 

=  (6t  +  24t  )  x  N  x  2cR 
m  a  -fy 


=  (6-  +  24t  )  x  (°1  °2  °3)  x  2cR 
2s  «■  ^ 


<6Tm  *  24Ta>  °1  °2  °3  2 
f64 


1-S 


c  R 


seconds 


(23) 


where  D^ ,  D2>  and  D3  are  the  maximum  spans  of  the  object,  and  s  =  0,  1,  or  2, 
as  discussed  in  the  previous  sub-section. 

From  Equation  23,  it  is  seen  that  the  central  processing  time  is  in¬ 
versely  proportional  to  the  frequency,  f,  and  to  the  fourth  power  of  the  lat¬ 
tice  cell  size,  6,  assuming  a  fixed  number,  R,  of  cycles  required  to  attain 
the  sinusoidal  steady  state.  The  major  implications  of  this  are  twofold: 


1.  Low-frequency  solutions  require  more  CP  time  than  high  fre¬ 
quency  solutions,  all  other  parameters  being  fixed;  and 

2.  Major  savings  of  CP  time  can  be  achieved  by  only  minor  increases 
of  6,  the  effective  spatial  resolution. 

Vector  Array  Processing  Machines 

Some  recent  computers  have  the  capability  to  perform  arithmetic  opera¬ 
tions  on  whole  arrays,  or  vectors,  of  numbers  in  a  highly  parallel  manner. 

Two  examples  of  such  machines  available  as  a  utility  service  are  the  Cyber 
203*  (Control  Data  Corporation  Cybernet)  and  the  Cray  1  (United  Computing 
Systems,  Inc.).  Such  a  capability  is  ideal  for  implementation  of  the  FD-TD 
method  since  virtually  all  of  its  operations  are  fully  explicit  and  are  per¬ 
formed  on  very  large  arrays  of  numbers. 

In  estimating  the  central  processing  time  for  the  FD-TD  method  on  an 
array  processing  computer  such  as  the  Cyber  203  or  the  Cray  1,  the  key  com¬ 
puter  parameter  is  the  machine  cycle  time,  x  ,  since  one  floating  point  multi¬ 
plication  or  addition  can  be  performed  by  such  machines  each  x  .  For  this 
case.  Equation  23  is  slightly  modified  to 

CP  »  30  xr  D,  D,  D,  21_S  cR  .  .  ,0/n 

c  1  2  3 _ seconds.  (24) 

f^4 

xc  is  approximately  12.5  nsec  for  the  machines  mentioned. 

6.3.2  Peripheral  Processing  Time 

At  each  time  step  of  the  FD-TD  algorithm,  the  following  twelve  data 
quantities  are  retrieved  from  computer  memory  for  each  unit  cell  of  the  space 
lattice:  1)  Three  vector  components  of  F;  2)  Three  vector  components  of  H; 
and  3)  Six  constant  coefficients.  Further,  the  following  six  data  quantities 
are  returned  to  computer  memory:  1)  Three  vector  components  of  F;  and  2) 
Three  vector  components  of  TT.  For  large  modeling  problems  involving  a  memory 
requirement  that  is  larger  than  the  maximum  central  memory  available  to  the 
user,  it  becomes  necessary  to  use  the  available  high-speed  mass  storage,  such 
as  magnetic  disc  files. 


♦successor  to  the  STAR-100 


In  this  case,  the  central  memory  is  used  only  as  a  scratchpad  or  buffer  to 
iteratively  access  portions  of  the  mars  memory,  perform  the  necessary  arith¬ 
metic  operations,  and  return  the  upd.  ced  electromagnetic  field  quantities  to 
the  mass  memory.  Data  transfer  can  be  programmed  either  in  an  explicit  fash¬ 
ion  or,  on  machines  with  virtual  memory,  can  be  performed  automatically  by 
the  computer. 

The  transfer  of  large  quantities  of  data  each  time  step  from  mass  stor¬ 
age  to  central  memory  and  back  to  mass  storage  requires  a  finite  length  of 
time  that  is  usually  placed  under  the  category  of  peripheral  processing  (PP) 
time.  Here,  the  key  factor  in  determining  PP  time  is  the  rate  of  data  trans¬ 
mission,  rrf  words/second,  to  and  from  the  mass  storage  device.  Using  the 
same  nomenclature  as  in  the  previous  subsections,  the  peripheral  processing 
time  for  one  time  step,  PPy  is  estimated  as 

PP,  a  18N  seconds  (25a) 

rd  time  step  ‘ 

The  peripheral  processing  time  for  the  modeling  of  the  entire  problem  is, 
from  Equations  25a  and  22b, 

PP  a  PP,  x  n 

1  max 

a  18N  x  2cR 
rd  f& 

a  18  x  (D1  °2  °3)  x  2cR 
rd  ^2S  6 3  ^ 

a  18  D1  P2  P3  21  S  cR  seconds.  (25b) 

rd  f«- 

Like  the  CP  time,  the  PP  time  is  inversely  proportional  to  the  frequency,  f, 
and  to  the  fourth. power  of  the  lattice  cell  size,  <5.  Again,  increasing  6 
only  slightly  helps  to  significantly  reduce  processing  time. 
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6.4  Impact  of  Requirements  for  Computer  Memory, 

CP  Time,  and  PP  Time _ 

During  the  present  research  effort,  many  large  FD-TD  programs  involving 
5  6 

N  >  10  cells  and  M  >  10  words  were  run  on  the  Control  Data  STAR-100  and 
Cyber  203  computers.  As  a  result  of  this  experience,  the  following  conclus¬ 
ions  have  been  reached  relative  to  the  impact  of  requirements  for  computer 
memory,  CP  time,  and  PP  time  in  large  FD-TD  models. 

1.  Limiting  M  to  less  than  the  central  memory  limit  of  the  computer, 
so  that  no  spill-over  onto  mass  storage  occurs,  results  in  rel¬ 
atively  very  small  amounts  of  PP  time.  Here,  the  total  program 
running  time  ("wall  clock  time")  is  essentially  equal  to  the  CP 
time,  which  is  in  the  range  of  a  few  minutes  for  the  Control  Data 
STAR- 100  or  Cyber  203. 

2.  If  any  mass  storage  is  used,  the  program  running  time  increases 
significantly  to  the  range  of  a  few  hours.  Here,  the  CP  time 
remains  very  low  for  the  STAR- 100  or  Cyber  203,  but  the  PP  time 
becomes  dominant.  Program  provisions  for  re-start  at  the  last 
completed  time  step  should  be  included  to  avoid  having  to  re-run 
the  entire  program  from  the  very  first  time  step  if  the  computer 
system  should  "go  down"  during  the  lengthy  program  running  period. 

3.  Usage  of  slower  machines  than  the  STAR-100  or  Cyber  203  would 
effect  mainly  the  CP  time.  Thus,  large  FD-TD  models  that  employ 
mass  storage  could  be  run  on  conventional  scalar  machines  such 

as  the  Cyber  76,  Uni  vac  1100,  or  Honeywell  6000  series  with  "wall 
clock"  times  similar  to  those  obtained  on  vector  processors  such 
as  the  STAR-100.  In  fact,  efficient  programming  of  the  accessing 
of  high  speed  disc  storage  devices  could  serve  to  reduce  "wall 
clock  time"  for  such  conventional  machines  and  provide  a  program 
re-start  capability. 

Overall,  in  determining  the  size  limits  of  an  object  that  can  be  modeled 
using  the  FD-TD  method,  one  must  consider  the  following  points: 
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1.  Will  mass  storage  be  used?  If  mass  storage  is  used,  there  is 
no  real  limit  to  M  of  Equation  21b,  so  that  arbitrarily  large 
objects  can  be  modeled,  in  principle.  However,  in  practice, 
the  object  size  is  limited  by  the  maximum  "wall  clock  time”  of 
the  program  which  is  acceptable  to  both  the  user  of  the  program 
and  the  operator  of  the  computing  center.  This  is  given  as  a 
function  of  the  spatial  dimensions  of  the  object,  the  lattice 
cell  size,  and  the  characteristics  of  the  computer  by  the  CP 
time  of  Equation  23  plus  the  PP  time  of  Equation  25b. 

2.  If  mass  storage  is  not  used,  the  value  of  M  of  Equation  21b  is 
limited  to  the  user-accessible  part  of  the  computer  central  mem¬ 
ory.  Now,  the  object  size  is  limited  sharply  by  Equation  21b 
because  larger  objects  (or  ones  with  higher  spatial  resolution) 
simply  cannot  be  fitted  into  the  computer  central  memory.  Here, 
the  question  of  increasing  the  lattice  cell  size,  6,  becomes 
important  because  of  the  strong  third-power  dependence  of  M  upon 
6. 
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7.0  PROGRAM  TASK  5:  INVESTIGATION  OF  INTERFACING  THE  FD-TD  METHOD 

_ WITH  OTHER  ANALYSIS  TECHNIQUES _ 

7.1  Introduction 

The  analysis  of  the  electromagnetic  excitation  of  an  aperture  on  an  arbi¬ 
trarily  shaped  object  is  a  complex  problem  [11].  This  problem  has  been  given 
special  attention  in  the  broad  area  of  EMP  penetration  and  simulation  studies 
[12],  with  great  efforts  being  applied  to  assess  coupling  to  objects  present 
behind  apertures  on  finite,  metallic,  hollow  scattering  bodies. 

The  problem  of  the  penetration  and  coupling  of  electromagnetic  energy 
through  an  aperture  has  been  studied  extensively  by  many  investigators  using 
analytical  and  method-of-moments  (MOM)  numerical  approaches  [13]  to  character¬ 
ize  the  behavior  of  simple  aperture  shapes  in  a  conducting  screen  [14]  -  [18] 
or  in  a  finite,  scattering  body  [19],  [20].  However,  the  analysis  becomes 
very  complex  if  there  are  other  scattering  objects  in  the  vicinity  which  are 
coupled  to  the  aperture  [11],  [21],  [22]  wherein  the  mutual  interaction  must 
be  fully  taken  into  account.  A  powerful  alternate  approach  is  the  pure  FD-TD 
method,  discussed  in  detail  in  the  previous  sections  of  this  report,  which 
allows  computation  of  the  penetrated  internal  electromagnetic  fields  by  direct 
modeling  of  the  object  of  interest,  including  apertures  and  internal  contents. 

In  order  to  treat  coupling  problems  more  effectively,  a  new  hybrid 
MOM/FD-TD  technique  has  been  developed  based  on  a  novel  use  of  the  so-called 
"field  equivalence  theorem"  due  to  Schelkunoff  [23].  This  hybrid  technique 
and  its  applications,  described  in  this  section,  is  basically  an  equivalent 
aperture  excitation  method.  This  allows  analysis  of  the  coupling  problem  in 
two  distinct  steps: 

Step  1  -  Analysis  of  the  relatively  simple  exterior  problem  using  the 
method  of  moments  (or  some  other  appropriate  technique)  to 
determine  the  equivalent  excitation  currents  and/or  fields 
in  the  aperture  regions.  This  can  be  done  independent  of 
any  knowledge  of  the  interior  contents  of  the  object  of 
interest. 
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Step  2  -  Use  of  the  FO-TD  method  to  analyze  the  relatively  complex 
interior  problem,  assuming  as  an  excitation  the  equivalent 
currents  or  fields  found  in  Step  1. 

In  this  way,  each  analysis  technique  is  applied  in  the  structure  size 
and  complexity  range  that  it  is  best  suited  for,  allowing  an  overall  solution 
that  is  accurate  for  large,  simple  structures  with  small,  complex  apertures 
or  other  fine  interior  details.  The  specific  details  of  adaptation  of  the 
data  obtained  from  the  MOM  frequency-domain  method  to  the  FD-TD  time-domain 
method  will  be  discussed  later  in  the  section  covering  hybrid  method  applica¬ 
tions. 

7.2  Basis  of  the  Hybrid  Technique 

In  this  sub-section,  the  hybrid  technique  based  on  the  use  of  the  field 
equivalence  theorem  by  Schelkunoff  [23]  (often  referred  to  as  Schelkunoff 's 
third  equivalence  theorem,  as  in  his  original  paper)  is  explained  in  a  sys¬ 
tematic  manner  to  analyze  aperture  field  penetration  and  coupling  into  struc¬ 
tures  by  interfacing  the  MOM  technique  (for  the  exterior-region  analysis) 
with  the  FD-TD  technique  (for  the  interior-region  analysis).  The  general  idea 
behind  the  hybrid  approach  is  not  new,  since  Schuman  and  Warren  applied  the 
Schelkunoff  third  theorem  to  interface  the  MOM  technique  with  MOM  to  separate¬ 
ly  analyze  the  exterior  and  interior  regions  of  an  open-ended  cylinder  and  a 
nose-cone  structure  with  a  nose  aperture  and  a  thin  circumferential  gap  aper¬ 
ture  [20].  Both  the  cylinder  and  nose-cone  structures  were  treated  as  bodies 
of  revolution.  Their  analysis  was  restricted  only  to  empty  interior  regions 
or  cavities.  If  the  interior  region  is  complicated  by  the  addition  of  inter¬ 
acting  material  bodies,  the  MOM/MOM  hybriding  becomes  very  complicated  and 
difficult.  For  such  a  situation,  a  MOM/FD-TD  hybriding  remains  possible,  as 
will  be  shown,  even  for  interior  material  bodies  strongly  coupled  to  the 
apertures. 

7.2.1  MOM  /  MOM  Hybrid  Technique 

Figure  33  illustrates  the  classical  problem  of  a  perfectly  conducting 
scatterer,  S^.,  with  an  aperture,  S^.  External  to  the  scatterer  in  Region  1, 
an  incident  field  (F1 ,  H1)  illuminates  the  aperture,  and  part  of  the  energy 
penetrates  into  the  cavity  (Region  2)  through  the  aperture.  To  compute  the 
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33.  PERFECTLY  CONDUCTING  SCATTERER  WITH  APERTURE 


81 


total  field  { ^ ,  Hg)  ’n  Re9i°n  2,  many  analytical  and  MOM  numerical  approaches 
studied  in  the  literature  have  employed  integral-equation  formulations  treat¬ 
ing  Regions  1  and  2  simultaneously.  These  approaches  characterize  either 
or  in  terms  of  tangential  fields  or  equivalent  surface  currents  [13]. 

Recently,  Schuman  and  Warren  [20]  have  developed  a  hybrid  technique  to 
treat  the  Region  1  problem  and  the  Region  2  problem  separately,  and  compute 
coupling  through  arbitrarily-shaped  apertures  in  conducting  bodies  of  revolu¬ 
tion.  In  this  manner,  rotational ly-symmetric  cavities  with  arbitrary  con¬ 
tours  can  be  modeled.  Further,  the  sensitivity  and  accuracy  of  field  compu¬ 
tation  is  greatly  improved  for  cases  when  the  aperture  coupling  is  small 
and/or  points  deep  within  the  cavity  are  sought.  The  latter  is  realized  be¬ 
cause  the  hybrid  method  avoids  any  substraction  of  nearly-equal  impressed 
and  scattered  fields  to  compute  the  total  field  within  Region  2.  Thus,  error 
due  to  residual  subtraction  noise  is  eliminated.  The  hybrid  approach  of 
Schuman  and  Warren  is  illustrated  in  Figure  34,  and  is  summarized  below  as  a 
five-step  computation. 

Region  1:  MOM  Technique 

(a)  The  aperture  region,  S^,  is  short-circuited,  and  the  straight¬ 
forward  exterior  problem  is  solved  [24],  [25]  to  obtain  the 

_j\ 

induced  electric  current  distribution,  J0,  in  the  short- 
circuited  aperture  region  (Figure  34b). 

(b)  The  short-circuit  current,  JQ,  is  now  placed  in  the  open  cir¬ 
cuited  aperture  region  with  a  sign  change  (Figure  34c)  to 
account  for  the  continuity  of  the  fields  in  the  aperture  region 
[23].  The  electric  current,  -J^,  fictitiously  placed  in  the 
open-circuited  aperture  region  acts  as  an  equivalent  source  for 
the  interior  fields  in  Region  2. 

Region  2:  MOM  Technique 

(c)  Based  on  the  generalized  multiport  approach  [13],  referring  to 

—A 

Region  2,  an  equivalent  electric  field  excitation,  E  ,  existing 
only  in  the  aperture  region  (Figure  34d)  is  determined  which 

A 

will  produce  the  electric  current,  -JQ,  with  the  aperture  short- 
circuited. 


REGION  (T)  ( §|  ,  H, ) 


(d) 

Fig.  34.  MOM/MOM  HYBRID  TECHNIQUE 
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(d)  The  electric  current  distribution,  J  ,  is  obtained  on  the  interior 

—A 

surface  of  the  conducting  body  with  E  as  the  source  equivalent 
electric  field  (Figure  34d). 

—  _  —C  —A 

(e)  The  interior  fields  (E0  H0)  are  now  computed  based  on  J  and  E  or 

_ A  £  9  <- 

■t- 

The  mathematical  formulation  of  the  hybrid  MOM  /  MOM  technique  [20], 

[24],  [25]  is  beyond  the  scope  of  this  report  and  is  not  treated  here.  How¬ 
ever,  one  point  should  be  repeated  from  the  original  paper  of  Schuman  and 
Warren  [20].  Namely,  any  loading  within  the  cavity  region  (Region  2)  changes 
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the  aperture  equivalent  electric  field  excitation ,  E  ,  of  Step  (c)  above. 
Therefore,  the  generalized  multiport  approach  of  Step  (c)  must  take  into  ac¬ 
count  any  loading  within  the  cavity  by  either  metal  or  dielectric  objects. 

This  requirement  introduces  considerable  complexity,  and  perhaps  explains 
why  the  MOM  /  MOM  hybrid  approach  has  not  been  more  widely  used  as  a  means  of 
modeling  realistic  problems.  This  cavity  loading  effect  is  discussed  in 
more  detail  in  Section  7.3,  where  it  is  shown  that  MOM  /  MOM  hybriding  is  very 
complicated  even  for  the  simple  case  of  aperture  loading  by  a  thin  wire  in 
the  cavity  region.  On  the  other  hand,  a  hybrid  formulation  employing  MOM  / 

A 

FD-TD  does  not  require  computation  of  E  ,  and  totally  avoids  the  complicating 
effect  of  loading  within  the  cavity. 


7.2.2  MOM  /  FD-TD  Hybrid  Technique 


In  order  to  determine  the  applicability  of  the  FD-TD  method  (a  time- 
domain  approach)  to  interface  effectively  with  other  numerical  procedures 
such  as  the  method  of  moments  (a  frequency-domain  approach),  the  hybrid  MOM  / 
FD-TD  technique  has  been  successfully  developed  and  implemented  to  compute 
aperture  coupling  and  field  penetration.  Again,  the  hybrid  formulation 
allows  one  to  treat  the  exterior  (Region  1)  problem  and  the  interior  (Region 
2)  problems  separately  in  two  distinct  steps.  In  Figure  35,  the  hybrid  MOM  / 
FD-TD  technique  is  illustrated  schematically,  and  is  summarized  below  as  a 
four-step  computation. 
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Fig.  35.  MOM  /  FD-TD  HYBRID  TECHNIQUE 


Region  1:  MOM  Technique 

(a)  Identical  to  Step  (a)  of  Section  7.2.1,  the  aperture  region, 

S^,  is  short-circuited,  and  the  straightforward  exterior  problem 

is  solved  to  obtain  the  induced  electric  current  distribution, 

Tg,  in  the  short-circuited  aperture  region  (Figure  35b).  In 

Appendix  A,  a  brief  summary  of  the  analytical  [24],  [25]  and 

—A 

numerical  formulation  [11],  [15],  [26]  needed  to  compute  JQ  on 
a  perfectly-conducting  body  of  revolution  is  discussed. 

(b)  Identical  to  Step  (b)  of  Section  7.2.1,  the  short-circuit  cur- 

i\ 

rent,  J0>  is  now  placed  in  the  open-circuited  aperture  region 
with  a  sign  change  (Figure  35c)  to  account  for  the  continuity  of 

_ A 

the  fields  in  the  aperture  region.  The  electric  current,  -0Q, 
fictitiously  placed  in  the  open-circuited  aperture  region  acts 
as  an  equivalent  source  for  the  interior  fields  in  Region  2. 

Region  2:  FD-TD  Technique 

(c)  We  note  that  the  MOM  technique  gives  the  equivalent  current 
source  distribution,  -Jq,  in  the  frequency  domain  as  a  spatial 
distribution  of  phasor  quantities  having  both  magnitude  (relative 

to  the  incident  fields)  and  phase  (relative  to  some  phase  reference, 
normally  at  the  origin  of  the  coordinate  system).  With  the  FD-TD 

method  being  a  time-domain  technique,  the  phase  distribution  of 

-A  .  ... 

-JQ  is  interpreted  as  a  time-delay  distribution  with  respect  to 

the  original  phase  reference  location.  The  magnitude  distribution 

is  taken  intact  without  interpretation  or  modification.  In  this 

manner,  the  FD-TD  aperture  equivalent  current  source  distribution 

assumes  sinusoidal  steady-state  quantities  starting  at  the  very 

beginning  of  time-stepping,  with  the  proper  time  delay  to  account 

for  phase  shift. 

(d)  Using  the  FD-TD  approach,  the  interior  fields  (£9,  H„)  are 

—A  ^  ^ 

computed  directly  by  using  -J0  as  a  source  term  distribution 
in  the  V  X  H  difference  equations.  No  additional  complication 
is  involved  if  the  cavity  is  loaded  in  any  manner,  since  this 
loading  is  modeled  directly  on  the  FD-TD  lattice.  Thus,  there 
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is  no  penalty  for  modeling  the  cavity  interior  to  as  fine  an 
extent  as  possible. 

7.3  Canonical  Case  Studies:  MOM  /  MOM  Hybrid  Technique 

In  this  section,  some  studies  are  reported  regarding  the  validity  of 
applying  Schelkunoff 's  third  equivalence  principle  for  apertures  having  wires 
closely  coupled  to  them.  Additional  results  are  also  reported  concerning  the 
need  to  model  the  entire  interaction  geometry  when  applying  the  equivalence 
principle,  instead  of  modeling  just  the  cavity  of  interest.  These  studies 
were  conducted  using  the  MOM  /  MOM  hybrid  technique  applied  to  the  canonical 
geometry  of  a  thin  wire  scatterer  placed  behind  a  narrow-slot  perforated 
conducting  screen,  as  shown  in  Figure  36.  Two  different  cases  were  examined: 

Case  A  -  Direct  straightforward  MOM  solution  of  the  entire  problem;  and 

Case  B  -  Hybrid  MOM  /  MCM  technique. 

These  will  now  be  discussed. 

Case  A:  Direct  straightforward  method-of -moments  solution 

The  general  theory  of  aperture  coupling  and  transient  characterization  of 
conducting  bodies  placed  behind  an  arbitrary  shaped  aperture  in  a  screen  has 
been  treated  in  Reference  11.  Specifically,  the  present  canonical  problem 
(Figure  36)  was  reported  in  Reference  27  with  numerical  results  derived  via 
the  direct  method  of  moments  solution.  Basically,  a  coupled  set  of  integro- 
differential  equations  [11],  [27],  [28]  was  derived  by  treating  the  aperture 
electric  field  distribution  and  the  wire  electric  current  distribution  as 
unknowns  by  enforcing  appropriate  boundary  conditions  at  the  aperture  and  on 
the  wire.  The  coupled  equations  were  solved  by  expanding  the  unknown  distrib¬ 
utions  in  terms  of  standard  pulse  functions  and  then  testing  the  resulting 
equations  using  piecewise  triangular  functions  [13].  A  brief  summary  of  the 
coupled  integral  equations  for  this  problem  is  given  in  Appendix  B. 

For  the  geometry  of  the  narrow  slot/thin  wire  shown  in  Figure  36,  we 
define  the  following  parameters: 
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Fig.  36.  GEOMETRY  OF  WIRE  SCATTERER  BEHIND  NARROW  SLOT 


SL=  length  of  narrow  slot  =  0.5  A 
w  =  width  of  narrow  slot  =  0.5  X 
L  =  length  of  wire  scatterer  =  0.5  X 
a  =  radius  of  wire  scatterer  =  0.001  X 
(xc>  y  ,  zc)  =  location  of  the  center  of  wire  scatterer  in  upper  half  space 
(z  >0) 

=  (0,  0,  0.25  X) 

6  =  angle  of  wire  relative  to  y  axis  =  0° 

Further,  the  narrow  slot  is  oriented  along  the  axis,  and  the  plane  wave 

->i 

(E  ,  H  )  excites  the  slot  region  in  the  lower  half  space  (z  <  0). 
y  x 

Figures  37  and  38  show  the  results  of  the  direct  method-of  moments 
solution.  Both  the  slot  and  the  wire  are  assumed  to  have  0.5  m  resonant 
lengths,  and  the  wire  is  placed  at  a  distance  0.25  m  directly  above  the  slot. 
The  slot  is  excited  by  a  300  MHz  plane  wave  in  the  lower  half  space  (z  <  0) 
at  normal  incidence  to  the  conducting  screen.  The  electric  field,  E^,  in  the 
slot  is  characterized  in  terms  of  an  equivalent  magnetic  current  distribution 


[27]  M  ,  in  the  shadow  side 

A 

A  A  A 

y  Ey  =  xMx  x  z 

(26) 

Mx  =  m(x)  £  (y) 

(27) 

where 


m(x)  =  axial  distribution  of  the  equivalent  magnetic  current  in  the  slot 

£(y)  =  transverse  distribution  of  the  equivalent  magnetic  current  in 
the  slot 

=  ^[(J)2  -  y2]_ls  (28) 

In  Figure  37,  the  shadow  side  (z  =  0+)  slot  axial  magnetic  current 
distribution  m(x)  in  the  presence  of  the  screen  and  the  wire  is  plotted  for 
6=0.  This  angle  corresponds  to  the  wire  axis  being  physically  perpendicular 
to  the  slot  axis,  which  is  the  orientation  for  maximum  coupling  and  inter¬ 
action  between  the  slot  and  wire.  Both  real  and  imaginary  parts  of  m(x) 
are  plotted.  Figure  38  plots  the  corresponding  distribution  of  the  induced 
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_ direct  method  of 

- moment  solution 

•  •  •  hybrid  solution 


w  =  0*05 
L  =0-5 
a  =0  001 
f  =300MI 


SHADOW-SIDE  SLOT  MAGNETIC  CURRENT  DISTRIBUTION 
IN  PRESENCE  OF  THE  SCREEN  AND  THE  WIRE 


Fig.  38.  WIRE  ELECTRIC  CURRENT  DISTRIBUTION  IN  THE 
PRESENCE  OF  THE  SCREEN  AND  THE  NARROW  SLOT 
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electric  current  on  the  wire.* 


Case  B:  Hybrid  MOM  /  MOM  technique 

The  limitation  of  Case  A  (which  involves  formulation  of  a  coupled  set  of 
integro-differential  equations  and  solution  via  the  straight  method  of 
moments)  is  that  one  has  to  analyze  the  complete  problem  in  a  simultaneous 
solution.  As  stated  in  Section  7.2.1,  an  alternative  is  to  decompose  the 
problem  according  to  the  MOM  /  MOM  hybrid  technique.  This  involves  the  suc¬ 
cessful  application  of  the  Schelkunoff  third  equivalence  theorem  [23]. 

To  apply  this  technique,  the  slot  in  the  canonical  geometry  shown  in 
Figure  36  is  first  short-circuited  by  replacing  the  slot  by  a  perfectly  con¬ 
ducting  material  identical  to  the  screen  itself.  Next,  the  lower  half  space, 
z  <  0,  is  analyzed  as  a  pure  boundary  value  problem.  The  induced  electric 

current,  J  ,  on  the  shorted  region  of  the  slot  is  determined  for  the  given 
o 

incident  excitation.  For  a  plane  wave  incident  normally  on  the  screen  with 

and  e""  oriented  along  x  and  -y  coordinate  axes,  the  equivalent  current 
x  y 

excitation  is  given  by 

J0  =  (29a) 

Next,  we  go  back  to  the  original  problem  and  place  in  the  open-circuited 
slotted  region  an  equivalent  aperture  electric  current  excitation  of 


=  2H’  y  .  (29b) 


*For  6  =  90°,  the  slot  and  wire  are  completely  uncoupled,  and  the  penetrating 
field  is  unaffected  by  the  presence  of  the  thin  wire.  Here,  the  m(x)  dis¬ 
tribution  has  been  verified,  since  it  is  the  Babinet  equivalent  of  a  narrow 
strip.  Further,  the  narrow  strip  distribution  can  be  verified  with  respect 
to  a  thin  wire  of  diameter  equal  to  half  the  strip  width. 

For  B  *  90°,  the  induced  electric  current  on  the  wire  is  zero,  representing 
complete  uncoupling  of  the  wire  and  slot. 
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The  sign  change  in  Equation  29b  is  in  accordance  with  the  boundary  condition 
requirements.  This  completes  the  first  stage  in  implementing  the  hybrid  tech¬ 
nique. 


The  second  stage  in  implementing  the  hybrid  technique  is  initiated  by 

-A 

computing  an  equivalent  aperture  electric  field  excitation,  E  ,  which  will 

produce  the  electric  current,  -J  ,  with  the  aperture  short-circuited.  When 

0  —A 

the  slot  and  wire  are  uncoupled  (8  =  90°),  E  is  found  by  treating  the  slot 
as  a  multiport  generalized  system  [13],  [20]  and  then  inverting  the  system 
equation : 

[3A]  =  [YA]  [EA]  (30a) 

where 

=A 

[Y  ]  =  slot  generalized  admittance  matrix. 


When  the  wire  is  oriented  such  that  8  1  90°,  there  is  significant  interaction 
between  the  wire  and  slot.  For  this  interesting  case,  the  multiport  general¬ 
ized  system  matrix  of  Equation  30a  must  be  modified  [13],  [20],  [28]  according 
to  Figure  39,  yielding 

[3A]  =  [[YA]  +  [YA]]  [EA]  (30b) 


where  [YA]  =  mutual  generalized  admittance  matrix  at  the  slot 
due  to  the  presence  of  the  wire  (accounts  for  wire 
coupling;  also  known  as  the  wire  admittance  loading 
matrix  or  transfer  admittance  matrix  [28]). 

and 

-A 

[Ew]  =  slot  equivalent  electric  field  excitation  in  the  presence 
of  the  wire  scatterer. 


For  the  simple  canonical  slot-wire  geometry  the  calculation  of  the 
[YA]  is  straightforward  [13]  and  in  fact  can  be  written  in  a  closed  form  for 
very  long  wire  scatterers  [28]. 
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Fig. 


39.  MULTIPORT  GENERALIZED  NETWORK  REPRESENTATION  FOR 
HYBRIDING  THE  METHOD  OF  MOMENTS  TECHNIQUE 
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The  hybrid  matrix  Equation  30b  has  been  solved  for  0  =  0°  in  the  canon¬ 

ical  geometry  of  Figure  36,  and  the  computed  results  for  the  slot  distribution 
are  plotted  in  Figure  37.  Excellent  agreement  with  the  direct  method-of- 
moments  solution  is  seen.  The  equivalent  aperture  electric  field  excitation 
radiates  in  the  presence  of  the  screen,  and  electric  currents  I  (c)  are  in¬ 
duced  on  the  wire  by  mutual  coupling  interaction.  In  Figure  38,  the  computed 
electric  current  induced  on  the  wire  (based  on  the  hybrid  formulation))is 
plotted.  Again,  excellent  agreement  with  the  results  obtained  for  the  direct 
method-of -moments  case  is  seen. 

Therefore,  it  is  concluded  that  it  is  valid  to  apply  Schel kunoff ' s  third 
equivalence  theorem  [23]  for  apertures  having  wires  closely  coupled  to  them. 

It  has  been  demonstrated  that  the  MOM  /  MOM  hybrid  approach  is  capable  of 
dealing  with  such  loaded  apertures  if  the  mutual  interaction  between  the  wire 
and  the  aperture  is  rigorously  taken  into  account.  The  next  section  will 
present  results  that  strongly  indicate  that  the  MOM  /  FD-TD  hybrid  technique 
is  also  valid  and,  in  fact,  even  much  more  useful  for  complex,  loaded  cav¬ 
ities  than  MOM  /  MOM  because  there  is  no  need  to  set  up  and  solve  for  mutual 
interactions  of  the  cavity  contents  and  the  apertures,  as  in  Equation  30b. 

7.4  Case  Studies:  MOM  /  FD-TD  Hybrid  Technique 

In  this  section,  three  case  studies  are  reported  regarding  the  validity 
of  the  MOM  /  FD-TD  hybrid  technique:  1)  Open-ended  cylinder,  axial  inci¬ 
dence;  2)  Open-ended  cylinder,  oblique  (6  =  45°,  <)>  =  0°)  incidence,  TM 
polarization;  and  3)  Loaded  missile  guidance  section,  axial  incidence.  The 
results  clearly  indicate  that  this  hybrid  technique  is  valid  and  useful  for 
these  and  similar  structures. 

7.4.1  Case  1:  Open-Ended  Cylinder,  Axial  Incidence 

The  MOM  /  FD-TD  hybrid  model  for  the  axial -incidence,  open-ended  aluminum 
cylinder  (19  cm  diameter  and  68.5  cm  long,  as  shown  in  Figures  40  and  41 )  was 
run  for  800  time  steps  (2.0  cycles  of  the  incident  wave  at  f  =  300  MHz).  The 
excitation  consisted  of  moment-method  data  derived  by  D.  Wilton  and  A.  Glisson 
([9]  and  Appendix  A)  for  the  electric  currents  over  the  aperture  locus  when 
shorted,  as  specified  by  Schel kunoff 's  third  theorem.  The  results  were 
contour-mapped  and  compared  to  previous  purely  FD-TD  results  obtained  in  RADC 
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Fig.  40.  VIEW  OF  MOM/FD-TD  HYBRID  CYLINDER  MODEL  AT  FD-TD 
LATTICE  PLANE  i  =  24.5  FOR  AXIAL  INCIDENCE  CASE 
(CASE  NO.  1) 
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Fig.  41.  CROSS-SECTION  OF  CYLINDER  MODEL  FOR  MOM/FD-TD 
HYBRID  CASE  NO.  1 
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Contract  F30602-77-C-0163. 

Figure  42  plots  the  comparison  of  the  pure  FD-TD  results  for  the  radial 

electric  field  along  the  cylinder  axis  with  the  pure  moment-method  results  of 

D.  Wilton  and  A.  Glisson  [9].  It  is  seen  that  the  results  agree  within  1  dB 
up  to  20  cm  back  from  the  aperture.  Beyond  this  point,  the  FD-TD  results  had 
not  quite  converged  after  800  time  steps. 

Figure  43  plots  the  comparison  of  the  pure  FD-TD  results  for  the  radial 

electric  field  along  the  cylinder  axis  with  the  hybrid  MOM  /  FD-TD  program. 

It  is  seen  that  the  results  agree  within  2  dB  up  to  40  cm  back  from  the 
aperture. 

Figures  44  to  49  compare  the  computed  field  contours  in  the  cylinder 
horizontal  and  vertical  symmetry  planes  (as  defined  by  Figure  41)  within  11 
cm  of  the  aperture.  Excellent  agreement  is  seen  for  all  contours  in  the  hor¬ 
izontal  plane,  and  for  the  Hx  contours  in  the  horizontal  plane.  Some  contour 
distortion  and  displacement  is  seen  for  the  Ez  and  E^  contours  of  the  vertical 
plane. 

Figure  50  and  51  plot  the  field  contours  computed  by  the  hybrid  program 
within  35  cm  of  the  aperture.  These  may  be  compared  to  Figures  18  and  17, 
respectively,  of  Report  RADC-TR-78-142  obtained  by  the  pure  FD-TD  approach. 

The  result  of  these  comparisons  is  the  conclusion  that  the  MOM  /  FD-TD 
hybrid  approach  yields  data  for  this  case  (a  0.19  X  circular  aperture  in  an 
elongated  cylinder)  that  are  consistent  with  previous  data  derived  from  both 
pure  MOM  and  pure  FD-TD  methods.  This  is  a  significant  agreement  of  three 
computational  techniques  based  upon  very  different  first  principles. 
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Fig.  45.  COMPARISON  OF  COMPUTED  Hx  CONTOURS  IN  HORIZONTAL 
SYMMETRY  PLANE 
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Fig.  47.  COMPARISON  OF  COMPUTED  Ez  CONTOURS  IN  VERTICAL 
SYMMETRY  PLANE 


104 


Distance  from  front 


105 


Distance  from  open  end  (cm) 


Distance  from  axis  (cm)  Distance  from  axis  (cm)  Distance  from  axis  (cm) 

(a)  Ez  (b)  Hx  (c)  Hy 

Fig.  50.  COMPUTED  FIELD  CONTOURS  IN  CYLINDER  HOR.  SYMMETRY  PLANE 
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Fig.  51.  COMPUTED  FIELD  CONTOURS  IN  CYLINDER  VERT.  SYMMETRY  PLANE 
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7.4.2  Case  2:  Open-Ended  Cylinder,  Oblique  Incidence, 

_ TM  Polarization _ 

The  MOM  /  FD-TD  hybrid  model  for  the  oblique-incidence  (45°),  TM  polar¬ 
ized  case  of  the  open-ended  aluminum  cylinder  (as  shown  in  Figures  52  and  53^ 
was  run  for  800  time  steps  (2.0  cycles  of  the  incident  wave  at  f  =  300  MHz). 
The  excitation  consisted  of  moment-method  data  derived  by  D.  Wilton  and 
A.  Glisson  ([9]  and  Appendix  a)  for  the  electric  currents  over  the  aperture 
locus  when  shorted,  as  specified  by  Schelkunoff 's  third  theorem.  The  results 
were  compared  to  pure  MOM  results  in  an  attempt  to  establish  validity  of  the 
MOM  /  FD-TD  hybrid  model  for  this  case. 

Figure  54  plots  the  comparison  of  the  hybrid  MOM  /  FD-TD  results  for  the 
axial  electric  field  within  the  cylinder  along  the  cylinder  axis  with  the  pure 
moment-method  results  of  D.  Wilton  and  A.  Glisson  [9].  It  is  seen  that  the 
results  agree  within  2.5  dB  in  the  span  3  cm  -  25  cm  back  of  the  aperture. 

This  is  nearly  the  same  level  of  agreement  observed  in  Figures  42  and  43  for 
the  axial -incidence  case.  This  implies  that  arbitrary  angles  of  wave  incid¬ 
ence  could  be  treated  using  the  hybrid  MOM  /  FD-TD  method  simply  by  specifying 
the  proper  aperture  equivalent  current  excitation,  without  having  to  change 
the  FD-TD  geometry  of  the  cylinder  (a  considerably  more  complicated  task). 
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Oblique  incidence,  TM  polarization  case: 
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Fig.  52.  VIEW  OF  MOM/FD-TD  HYBRID  CYLINDER  MODEL  AT  FD-TD 
LATTICE  PLANE  i  =  24.5  FOR  OBLIQUE  INCIDENCE, 
TM  POLARIZATION  CASE  (CASE  NO.  2) 
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Fig.  53.  CROSS-SECTION  OF  CYLINDER  MODEL  FOR  MOM/FD-TD 
HYBRID  CASE  NO.  2 
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Fig.  54. 


COMPARISON  OF  RESULTS  FOR  THE  AXIAL  ELECTRIC  FIELD 
WITHIN  THE  CYLINDER,  A5°/TM  INCIDENCE  CASE 


7.4.3  Case  3:  Loaded  Missile  Guidance  Section, 

_ Ax  i  a  1 _ Incidence _ _  _ 

The  MOM  /  FD-TD  hybrid  model  for  the  axial-incidence  case  of  the  loaded 
guidance  section  (as  shown  in  Figures  55  and  24-26)  was  run  for  1800  time 
steps  (3.0  cycles  of  the  incident  wave  at  f  =  300  MHz).  The  excitation  con¬ 
sisted  of  moment-met  hoc  data  derived  by  D.  Wilton  and  A.  Glisson  ([9]  and 
Appendix  A)  for  the  electric  currents  over  the  aperture  and  sleeve-fitting 
loci  when  shorted,  as  specified  by  Schel kunoff ' s  third  theorem.  The  results 
were  compared  to  pure  FD-TD  results  presented  in  Section  4.2  of  this  report 
in  an  attempt  to  establish  the  consistency  of  the  MOM  /  FD-TD  hybrid  model 
for  this  case. 

Figure  56  plots  the  comparison  of  the  hybrid  MOM-FD-TD  results  for  the 
H  field  contours  in  the  vertical  symmetry  plane  with  the  pure  FD-TD  results 
already  presented  in  Figure  28.  It  is  seen  that,  for  botn  methods  stepped  to 
1800  time  steps,  there  is  an  excellent  agreement  of  the  0  dB  contours  near 
the  wire  connecting  the  cooled  detector  unit  to  the  pre-amp  can.  For  this 
contour,  the  maximum  spatial  shift  is  only  about  0.16  cm  in  a  direction 
further  out  from  the  wire.  Further,  there  is  observed  to  be  excellent  agree¬ 
ment  of  the  +10  dB  contour  near  the  wire  connecting  the  pre-amp  can  to  the 
metal  backplane.  For  this  contour,  the  maximum  spatial  shift  is  only  about 
0.1  cm  in  a  direction  further  out  from  the  wire.  This  implies  that  the  cur¬ 
rents  in  these  two  major  wires  are  predicted  to  be  almost  the  same  by  both 
the  pure  FD-TD  method  and  the  hybrid  MOM  /  FD-TD  method,  following  the  dis¬ 
cussion  of  Section  5  of  this  report. 

Figures  57  and  58  plot  the  comparison  of  the  hybrid  /  FD-TD  results 
for  the  Hx  and  Ez  fields  along  a  vertical  cut  through  the  center  of  the 
guidance  section  at  a  point  21  cm  in  back  of  the  nose  aperture  (about  2  cm  in 
front  of  the  sleeve  fitting,  at  the  point  where  the  circumferential  slot 
opens  into  the  interior  of  the  nose  cone).  The  hybrid  run  results  are  after 
1500  time  steps,  while  the  pure  FD-TD  results  are  after  1800  time  steps.  For 
this  case,  a  very  high  level  of  agreement  is  observed  between  the  two  sets...cf 
data  at  all  points  of  comparison.  The  worst-case  difference  between  the 
results  is  only  1  dB,  with  most  results  consistent  within  only  fractions  of 


f  =  300  MHz 


Fig.  55.  VIEW  OF  MOM/FD-TD  MISSILE  GUIDANCE  SEC! ION  MODEL 
AT  FD-TD  LATTICE  PLANE  i  =  24.5  FOR  AXIAL 
INCIDENCE  (CASE  NO.  3) 
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of  a  decibel.  (Comparison  of  the  data  sets  for  the  1800  time-step  case  for 
each  set  results  in  slightly  lessened  agreement,  such  that  the  worst-case 
difference  is  about  1.7  d8.  This  may  result  from  the  hybrid  program  progres¬ 
sing  to  the  sinusoidal  steady  state  at  a  slightly  faster  rate  than  the  pure 
FD-TD  program,  since  a  sinusoidal  steady-state  equivalent  aperture  current 
excitation  is  employed  from  the  very  beginning  of  the  hybrid  program,  rather 
than  an  aperture  excitation  which  must  build  to  the  steady  state.) 

This  case  study  leads  to  the  following  observations  and  conclusions: 

1.  The  hybrid  MOM  /  FD-TD  approach  yields  results  for  the 
missile  guidance  section  wire  currents  and  electromagnetic 
fields  which  are  consistent  with  the  pure  FD-TD  data. 

2.  For  the  axial -incidence  case,  the  assumption  in  Section  4  of 
this  report  that  the  finite  length  of  the  missile  has  little 
effect  on  the  guidance  section  fields  is  shown  to  be  valid, 
since  the  hybrid  case  aperture  excitation  does  take  into 
account  the  missile  length,  and  its  results  are  consistent 
with  those  of  Section  4. 

3.  A  wire  passing  very  close  to  an  aperture,  and  strongly  coupled 

to  that  aperture  (as  for  the  case  of  the  pre-amp  can-to-backplane 
wire  near  the  sleeve-fitting  aperture)  can  be  consistently 
modeled  using  the  pure  FD-TD  approach  and  the  hybrid  MOM  /  FD-TD 
approach.  In  the  latter,  it  must  be  remembered  that  the  M0M- 
derived  equivalent  aperture  excitation  takes  into  account  none 
of  the  interior  details  of  the  structure. 

4.  Complex  sub-sections  of  elongated  simply-shaped  structures  are 
prime  candidates  for  detailed  modeling  of  the  interior  penetrating 
fields  via  the  hybrid  MOM  /  FD-TD  method. 

MOM-derived  sinusoidal  steady-state  equivalent  aperture  excitations 
can  be  easily  interpreted  to  provide  input  data  for  the  FD-TD  method. 
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8.0  DISCUSSION  AND  CONCLUSIONS 


This  research  program  demonstrated  that  the  FD-TD  method  can  be  success¬ 
fully  applied  to  electromagnetic  coupling  problems  involving  highly-conducting 
structures  with  hole  and  sleeve-type  apertures  illuminated  by  a  plane  wave 
having  an  arbitrary  polarization  and  angle  of  incidence.  Further,  the  FD-TD 
method  can  be  applied  to  complex  cavity-like  structures  having  internal  metal 
and  dielectric  materials,  as  well  as  connecting  wires.  Accuracy  of  the  FD-TD 
results  was  very  good  relative  to  the  uncertainties  of  avilable  experimental 
and  numerical -theory  approaches.  Convergence  of  the  electromagnetic  fields 
to  the  sinusoidal  steady  state  occurred  within  about  3  cycles  of  the  incident 
wave  when  a  slight  value  of  isotropic  conductivity  was  assigned  to  the  inter¬ 
ior  of  the  structures  modeled.  This  resulted  in  program  central  processor 
times  of  less  than  5  minutes  for  FD-TD  lattices  containing  as  many  as 

1.5  million  unknowns  time-stepped  to  nm„„  =  800  using  the  Control  Data 

r  max  3 

STAR-100  and  Cyber  203  computers. 

This  research  program  also  established  the  feasibility  of  a  hybrid 
MOM  /  FD-TD  analysis  technique  based  upon  the  use  of  Schel kunoff ' s  equiv¬ 
alent  current  theorem  at  apertures  of  a  cavity.  This  hybrid  technique  was 
shown  to  give  consistent  results  for  structures  illuminated  at  both  oblique 
and  axial  incidence  by  plane  waves.  Further,  this  technique  was  shown  to 
give  consistent  results  for  apertures  strongly  coupled  to  internal  wires  or 
other  metal  and  dielectric  structures. 

The  pure  FD-TD  and  hybrid  MOM  /  FD-TD  methods  appear  to  have  great 
promise  for  applications  involving  complex  conducting  and  dielectric  struc¬ 
tures  illuminated  by  plane  waves  at  arbitrary  angles  of  incidence  and  polar¬ 
ization.  Further,  great  promise  is  shown  for  those  structures  that  are  sim¬ 
ultaneously  electrically  large  compared  to  a  wavelength  and  penetrated  by 
locally  complex  cavity-backed  apertures  having  dimensions  comparable  to  a 
wavelength.  Such  structures  requiring  resolution  of  both  large  and  small 
details  at  the  same  time  have  not  been  well  treated  by  any  one  previous 
analytical  or  numerical  approach. 
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APPENDIX  A 


Analytical  and  Numerical  Formulation  of  the 
Perfectly  Conducting  Body  of  Revolution 


This  appendix  summarizes  the  analytical  formulation  and  the  MOM  num¬ 
erical  procedure  to  treat  scattering  by  perfectly  conducting  bodies  of  revo¬ 
lution  [24]  -  [26].  Due  to  the  symmetry  of  such  structures,  the  incident 
plane-wave  excitation  and  corresponding  induced  electric  current  distribution 
are  expanded  in  terms  of  a  sum  of  Fourier  modes  to  reduce  the  coupled  integro- 
differential  equation  to  a  simpler  modal  generalized  matrix.  This  matrix 
can  then  be  solved  for  every  mode  excited  on  the  body  of  revolution. 

I.  Integral  Equation:  Coupled  Form 

The  geometry  of  a  perfectly  conducting  body  of  revolution  is  shown  in 
Figure  A-l .  It  is  placed  in  a  homogeneous  free  space  medium  of  permeability, 
UQ,  permittivity,  eQ,  and  conductivity,  aQ  =  0.  The  z  coordinate  axis  of  the 
(p,  <j>,  z)  coordinate  system  coincides  with  the  axis  of  the  body  of  revolution. 
P1  is  a  point  (p‘,  <p 1 ,  z‘)  on  the  surface  of  the  body  of  revolution  with  t' 

A  A 

and  n'  unit  tangential  and  normal  vectors.  The  unit  tangential  vector,  t‘, 

AAA 

makes  an  angle,  y' ,  with  the  z  axis.  Throughout  the  analysis,  p  t  x  n 
(right  hand  system)  is  assumed. 


[En(r),  H^r)]  =  external  incident  plane  wave 

J(r)  =  induced  surface  electric  current 

^  distribution 

[Es(r),  Hs{r)]  =  scattered  electric  field 


On  the  surface  of  the  body  of  revolution 
condition  is  satisfied, 


n  x 


EV)  = 


-  n  x 


ES(n) , 


r  c  on  S 


S,  the  following  boundary 

(A-l) 


Further,  the  scattered  electric  and  magnetic  fields  are  given  by 

ES(r)  =  -  j<a  A(r)  -  V  $(r)  (A-2a) 

Hs(r)  =  1  7  x  A(r)  (A-2b) 
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where  the  magnetic  vector  potential  ^(r)  and  the  scalar  potential  <j>(r)  can  be 
expressed  in  terms  of  the  surface  electric  currents  and  charges  on  the  body 
of  revolution: 


A(r) 


//;<?■, 


G  (r,  r1)  ds ' 


(A-3a) 


(A-3b) 


where  the  free-space  Green's  function  is 


G(r,  ?')  =  e‘jkR 

R 

R  =  [p2  +  p'2  -  2  pp*  Cos  (4>  -  4>')  +  (z  -  z')2]1* 


(A-3c) 

(A-3d) 

(A-3e) 


In  Equation  A-3c,  k  is  the  propagation  constant  in  the  free  space  medium,  and 
a)  is  the  frequency  of  the  incident  excitation  in  radians/sec.  In  obtaining 
Equations  A-2  and  A-3,  the  field  quantities  are  assumed  to  have  an  eJwt  time 
dependence.  Hence,  on  substituting  Equations  A-2  and  A-3  into  the  boundary 
condition  of  Equation  A-l ,  the  following  integro-differential  equation  is 
obtained: 

n  x  E  (r)  =  n  x 


d(r' )  G(r,  r' )  ds ' 


4ttojc, 


U 


Pe(r')  G(f,  r’) 


ds'  , 


S 


r  c  on  S 


(A-4) 


in  which  we  can  replace  the  charge  distribution  by 

p6(7,)  “  j  [v;.J(r')] 


(A-5) 
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In  Equation  A-4,  the  unknown  electric  current  distribution  J(r),  can  be 

A 

split  into  two  orthogonal  components  along  t  and  <J>  to  utilize  the  symmetry  of 
the  structure  in  the  <j>  coordinate  variable.  Similarly,  the  incident  field 
can  be  resolved  into  two  orthogonal  components.  This  gives: 


n  x  J(r)  =  t  Jt(r)  +  0  J^r)  (A-6a) 

n  x  E^r)  =  t  E^(r)  +  $  E^(r)  (A-6b) 

where  the  unit  vectors  are  given  by 

A  A  A  A 

t  =  sin  y  cos  <{>  x  +  sin  y  sin  $  y  +  cos  y  z  (A-7a) 

A  A  -A  /A 

n  =  cos  y  cos  4>  x  +  cos  y  sin  <p  y  -  sin  y  z  ( A- 7b ) 

-A  A  A 

<|>  =  -sin  <|>  x  +  cos  <p  y  (A-7c) 


Hence,  separating  the  t  and  <j>  components  of  Equation  A-4  gives  the  following 
two  coupled  integral  equations: 


M2 

(A-8a) 

L22  <V 

(A-8b) 

r  e  on  S 


where Ljnn  is  an  operator  given  by  the  integral  expressions 

Ln  (Jt>  =  ^ '0 JJJt(?,)  G(r,  r')  ds1 

S 

S 

L)2<v  =  ^JJv?,>  G< 


r ,  r ' )  ds' 


+  J  2 

4ttw£  9t 
0 


m 


=3 


[J  fr1)]  G(r,  r’)  ds1 


(A-8c) 


(A-8d) 
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L21  (Jt}  =  J^° jj Jt (;,)  c^,]  G(;* ?,)  ds' 

s 

fffr  Jt(?')]]  «(?.?■ 


4™eQp  d<t> 


)  ds’ 


M  \  =  JWU 

'22  vV  ' 


J^(r' )  Cd>*4>‘ ]  G(r,  r' )  ds1 


&J7 

S 

iff  ft-  ir[Vf-?'^G(f-f,)ds' 


where 

/S  A. 

ft  =  sin  y  sin  y'  cos  (<t>'~<p)  +  cos  y  cos  y' 

/\  a 

ft  -  sin  y  sin 

A  A 

4> - 1 ’  =  sin  y'  sin  (c^ 1  -cp ) 

A  A 

cj)  •  <p 1  =  cos  (<t>‘  -<p) 


(A-8e) 


(A-8f ) 


(A-8g) 
(A-8h) 
(A-8i ) 
(A-8j) 


II .  Modal  Form  of  Equations 


The  cp  symmetry  of  the  body  of  revolution  can  be  utilized  to  reduce  the 
double  integral  of  Equations  A~8c  to  A-8f  to  single  integrals  by  performing 
first  the  integration  with  respect  to  the  <)>'  variable.  For  this,  the  kernel 
function  G(r,  r' )  is  expanded  in  a  Fourier  series  as  follows: 


G(<p-cp 1 )  =  1  y  r  jm( <p-<p ‘ ) 
2?  m  e 

m=-co 

where 

Gm  =  G-m  =  /  GU)  cos  (mC)  dC 

-tt 


(A-9a) 


(A-9b) 
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in  which  G(£)  is  evaluated  on  the  surface  S, 

-JkR0 
G(C)  =  e _ ° 


K  =  [p2  +  p'2  -  2pp'  cos  e  +  (z-z')2]h 


Further,  Equation  A-9a  yields 


(A-9c) 


(A-9d) 


cos  ' )  S(4>-4> ' ) 

sin  ( 1 )  G((j>-<j)' ) 


where 


G 


c 

m 


rGm-l 

{rZT 


6m+l 
2j  J 


=  E 
m=-°° 


=  E 
m=-°° 


Gc  ejm(<t>-4>‘ ) 
m 

GS  ) 


(A-9e) 

(A-9f ) 

(A-9g) 

(A-9h) 


The  Fourier  expansions  of  Equations  A-9a,  A-9e,  and  A-9f  are  now  sub¬ 
stituted  into  the  integrals  of  Equations  A-8a  and  A-8b.  These  equations  are 
then  tested  on  both  sides  with  respect  to  the  ^-coordinate  variable  on  the 
body  of  revolution  by  multiplying  throughout  by  e"JP<*>  and  integrating  between 
the  limits  -tt  to  tt.  The  following  integral  equations  are  obtained,  with  the 
orthogonal ity  relationship  that  the  integrals  exist  only  for  m  =  p: 

/  Ej  (t,  4>)  e-jm*  d4>  (Jt)  +L12  (A-10a) 

-TT 

/  Ej(t,  ♦)  d«  -L2,  (Jt>  +U2  <V 

“IT 
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where 


Ln  fJt)  -  J“u0  J  J  Jt  (f.  *')  Dm  (p'dt’dd1 ) 


trA_X  ff 

47 uoe^  3t  J  J 
0  f  4>' 


3  xi  m  ti  a~  3m<f) 


4)')]  Gme-jnKp  (dfdd)1) 


(A-lOc) 


'V  ■  £ 


~  // °*  (t’’  E"  (p'dt'd*') 

t'  <f>' 

hr  h  ff  w  V"'(dt'dd') 


47Toje0  3t 


(A-1CJ) 


L*  (Jt’  ^  "■’ 


Fm  e‘jm<!!  (p'dt'd*') 


t'  <|)' 


+  4iir  f  f  lkr  [p‘  jt(t'  *  )]  (jm)  Gme  j  ^  (dt‘d*' ] 
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where 


D  =  sin  y  sin  y'  G*"  +  cos  y  cos  y'  G 


E  =  sin  y  Gm 
m  m 


F  =-sin  y'  Gm 
m  m 


(A-lOg) 


(A-lOh) 
(A-lOi ) 


(A-lOj) 
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m  m 


The  surface  electric  current  distribution  directed  along  t  and  <p  can 
also  be  expanded  in  terms  of  Fourier  modes,  as  follows: 


Jt  (t’’  ♦’>  =  2^-Z  Jt  «">  6 


jn<p' 


n=-oo 


(A-lla) 

(A-llb) 


n=-°° 


Since  the  Gm  coefficients  in  Equation  A- 1 0  depend  only  on  (t,  t‘,  z,  z'),  the 
above  current  expansions  can  be  substituted  into  the  integrals  of  Equations 
A-lOc  to  A-lOf.  Again  invoking  the  orthogonality  relationship  that  the  inte¬ 
grals  exist  only  for  m=n  with  respect  to  the  tp'  integration,  we  further  ob¬ 
tain  the  following  coupled  integral  equations  in  terms  of  only  ( t * )  and 
( t ' )  along  a  particular  generating  arc  of  the  body  of  revolution: 
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(A- 1 2d ) 
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III.  Reduction  to  Modal  Matrix  Form 

The  smooth  generating  arc  is  now  approximated  by  small  linear  segments 
as  shown  in  Figure  A-2.  In  Equations  A- 12a  and  A-12b,  the  unknown  currents, 
J™  ( t ’ )  and  ( t ' ) ,  are  expanded  by  piecewise  pulse  functions  shown  in 
Figures  A-2a  and  A-2b: 

N 

J?  (f  )  =  T  I?”  P?  (f)  ( A-  13a) 

1  z  1 


P?  (f)  =  i1  tn  t  <  f  <  t  ( A- 1 3b ) 

i  q  n-4  ~  ~  "+% 

(t’)  =  T  l|IA  p}  (t1)  (A-13c) 

*  7*1  ♦  2 
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we  obtain  the  linear  algebraic  modal  equations: 
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n  =  1,  2,  . N 

v  =  1 ,  2,  . N 

t  =  1 ,  2,  . N+l 

p  =  1  ,  2,  . N+l 

where  the  various  modal  impedance  parameters  are  given  by, 
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and 
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Now,  the  left-hand-side  excitation  terms  in  Equations  A-15a  and  A-1 5b  are 
obtained  for  a  particular  incident  plane-wave  excitation.  For  negative 
values  of  the  modes  (-m),  Z,-,  =  Z,,  Z,7  =  -Z,7  Z71  =  -Z01  and 

I  lv  I  lv,  I4V  I4V,  i.  ip  t  ip 

7-m  Z  7mZ 

LZZn  "  l22‘ 

P  P 
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IV.  Modal  Excitation  Terms 


In  Equation  A-4,  for  a  plane  wave  incident  on  a  body  of  revolution 
(Figure  A-l ) ,  we  have 


..  -*■ 
_  ri  Q-Jk  n  -r. 


E1  (r)  =  e 

and  in  terms  of  6-  and  <j>-  polarization  of  the  incident  field  we  have 


(A-16a) 


Ei  ■  Ee  6'  +  EI  ♦’  ■ 

The  plane  wave  can  be  resolved  into  a  sum  of  cylindrical  modes  given  by 

g-jkn1 . r  =  ejkz  cos  e1^-  sin  Qi }  ejn(4>V )  (A-16c) 

n=-oo 

where  J  is  the  Bessel  function  of  order  n.  According  to  Equations  A-l 2a  and 

^  i  i 

A-12b,  by  performing  the  modal  expansion  for  the  Et  and  E^  components  of  the 

plane  wave  excitation,  we  obtain  the  following  driving  terms  for  Equations 

A-15a  and  A-l 5b: 
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and 


VJp  =  [Ee  cos  01  jm  CJm-1(kpp-i5  sin  01 )  +  Jrn+1  (kpp-i5  sin  01 )] 
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where  Eg  and  are  the  amplitudes  for  the  0-  and  <p-  polarizations.  We  note 
in  Equations  A-17a  and  A- 17b  the  driving  vector  for  a  given  mode  is  the  sum 
of  contributions  from  the  two  0-  and  4>-  polarizations.  For  negative  values 
of  the  modes  (-m),  Vi"10  =  vf ,  V"m0  =  -V™0,  and  v;m<p  =  -V?*,  VT"1*  =  V?. 

t,.  t,.  <p„  <P_,  t  t.  <J>„  <J>_ 


v  Tp 


4 

_v  Tp 


V.  Modal  Matrix  Solution 

According  to  Equations  A-15a  and  A- 15b,  we  have  the  matrix  equation  for 
mode  m, 
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The  matrix  equation  is  solved  and  the  solution  for  the  induced  current 
vectors  obtained  both  for  0-  and  4>-  polarizations,  for  +m  modes 
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and  for  -m  modes 
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Hence,  referring  to  Equations  A-lla  and  A-llb;  A-13a  and  A-13b;  and  A-19a  and 
A-19b,  the  distribution  of  the  induced  current  for  a  given  plane  wave  excita¬ 
tion  is  given  by: 


1  N  °° 
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VI.  Application  to  the  Hybrid  MOM/FD-TD  Case  Studies 
of  Section  7.4 _ _ _ 

Cases  1  and  2:  Illumination  of  an  Open-Ended  Cylinder 

Figure  A-3  depicts  the  coordinate  systems  used  for  the  MOM/FD-TD  hybrid 
analysis  of  the  open-ended  cylinder  for  both  Case  1  -  Axial  Incidence  and 
Case  2  -  Oblique  (45°)  Incidence  with  TM  polarization.  Note  the  different 
coordinate  systems  used  for  the  MOM  part  of  the  analysis  [9]  and  the  FD-TD 
part  of  the  analysis. 

Figure  A-4  graphs  the  magnitude  of  the  MOM  modal  coefficients*  of  the 
equivalent  surface  electric  current  at  the  cylinder  aperture  for  Case  1. 
Since  only  the  m=l  mode  is  present,  the  equivalent  aperture  current  is  given 
by  [9] 

A  A 

0  (t,  <(>)  =  t  (t)  cos  <j)  +  4)  j  (t)  sin  <p  amps/meter  (A-21 ) 
for  any  point  (t,  <t>)  on  the  aperture.  This  can  be  easily  converted  to  the 

A  /V 

x  and  z  coordinates  of  the  FD-TD  system  to  derive  Jx  and  Jz  for  step  2  of 
the  analysis,  as  follows: 

A  A 

Jx  =  (4>  •  J)  cos  <j>  -  (t  •  J)  sin  <J>  amps/meter  (A-22a) 

A  /V 

0  =  (<|>  •  J)  sin  <J>  +  (t  •  0)  cos  $  amps/meter  (A-22b) 

At  all  locations  of  E  and  E  in  the  cylinder  aperture,  J  and  J  can 

A  4  A  L. 

now  be  inserted  as  field  sources  via  the  V  x  H  part  of  the  FD-TD  algorithm. 
Specifically,  Equations  9d  and  9f  of  Section  2.2.1  are  re-written  for  these 
points  as 


*For  clarity,  only  the  magnitude  is  graphed.  However,  in  computing  J  ( t ,0 ) , 
both  the  real  and  imaginary  parts  are  required. 
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FD-TD  z  axis 


t,  <t>:  method-of-moments  program  co-ordinates 

J:  surface  current  derived  from  method-of-moments  program 


Fig.  A-3 


COORDINATE  SYSTEMS  FOR  MOM/FD-TD  HYBRID  ANALYSIS 
OF  OPEN-ENDED  CYLINDER  AT  CYLINDER  APERTURE 
(AXIAL  INCIDENCE  CASE  AND  45° /TM  INCIDENCE  CASE) 
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-4.  MODE  COEFFICIENTS  FOR  APERTURE  EQUIVALENT  ELECTRIC 
CURRENT,  HYBRID  MOM/FD-TD  FORMULATION  OF  AXIAL 
INCIDENCE  ILLUMINATION  OF  OPEN-ENDED  CYLINDER 
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m  =  MEDIA(i+*s,j,k) 


^n+1  _n  n+h  n+Ji 

Ex(i+is,j,k)  =  ca(m)Ex(i+i5,j  ,k)  +  Cb(m)[Hz(i+isJ+^,k)  -  Hz(i+4,j-4,k)  + 

n+h  n+4  n+% 

Hy(i+ii,j  ,k-4)  -  Hy(i+i$,j  ,k+4)  +  v3 x ( i  +i^ , j  , k ) ]  (A-23a) 

m  =  MEDIA(i,j,k+%) 

~n+l  n+Jj  n+% 

Ez(i,j,k+%)  =  Ca(m)Ez(i ,j ,k+4)  +  cb(m)[Hy(i+Js,3 Mh)  -  H  (i-4,j ,k+4)  + 

n+k  n+*s  n+J$ 

H  (i.j-^.k+'-s)  -•  H  (i  ,j+4,k+4)  +  J  (i  ,j  ,k+H)]  (A-23b) 

A  X  Z 

In  Equation  A-23a,  specifying  Jx  in  the  aperture  is  equivalent  to  spec¬ 
ifying  the  addition  of  a  discontinuity  in  the  z-directed  tangential  magnetic 
field,  Hz,  across  the  aperture  source  plane.  Namely,  we  have  added 

Jx  =  H'  -  Hz,  (A-24a) 

where  H*  and  H*  are  tangential  magnetic  fields  located  an  infinitesimal  dis¬ 
tance  to  either  side  of  the  equivalent  aperture  source  plane.  This  procedure 
is  consistent  with  that  originally  discussed  by  Schelkunoff  [23].  In  a  sim¬ 
ilar  manner,  specifying  Jz  in  the  aperture  in  Equation  A-23b  is  equivalent  to 
specifying  the  addition  of  a  discontinuity  in  the  x-directed  tangential  mag¬ 
netic  field,  Hx,  across  the  aperture  source  plane,  namely 

J2  =  H“  -  Hx.  (A-24b) 
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Figures  A-5a  and  A-5b  graph  the  magnitude  of  the  MOM  modal  coefficients 
of  the  equivalent  electric  current  at  the  cylinder  aperture  for  Case  2.  Here, 
it  is  seen  that  modes  0,  1,  2,  and  3  predominate  for  J^.,  and  modes  1,2,  and 
3  predominate  for  Thus,  the  equivalent  aperture  current  is  given  by  [9] 


t  =  t  {Jf  +  E  J.  cos  (m<j>)} 

1  Lo  m=l  zm 

3 

=  -p  {Jt  +  E  cos  (nr»4>) }  (A-25a) 

o  m=l  ^m 


J.  <C  =  <t>  {E  J  sin  (m4>)} 

*  m=l  % 

.  3 

-  <t>  (E  0.  sin  (md> ) >  (A-25b) 

m=l  9m 

A  A 

for  any  point  (t,  <p)  on  the  aperture.  (Note  that  t  -  -p  on  the  aperture 

A  A 

locus.)  Again,  this  can  easily  be  converted  to  the  x  and  z  coordinates  of 
the  FD-TD  system  to  derive  Jx  and  0z  for  step  2  of  the  analysis  using 
Equations  A-22a  and  A-22b. 
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Fig . A-5a .  MODE  COEFFICIENTS  FOR  APERTURE  EQUIVALENT  ELECTRIC 
CURRENT,  HYBRID  MOM/FD-TD  FORMULATION  OF  OBLIQUE 
INCIDENCE  TM  ILLUMINATION  OF  OPEN-ENDED  CYLINDER 
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Case  3;  Illumination  of  a  Missile  Guidance  Section 

Figure  A-6  depicts  the  geometry  and  coordinate  systems  used  for  the  MOM/ 
FD-7D  hybrid  analysis  of  the  loaded  missile  guidance  section.  The  coordinate 
systems  used  are  the  same  as  shown  in  Figure  A-3  for  the  open-ended  cylinder. 

Figure  A- 7  graphs  the  magnitude  of  the  MOM  modal  coefficients  of  the 
equivalent  electric  current  at  the  nose  cone  sleeve-fitting  aperture  and  the 
nose  aperture.  Similar  to  the  axial-incidence  case  of  the  cylinder,  only  the 
m»l  mode  is  present. 

For  the  nose  aperture,  both  and  coefficients  were  used  to  form 
J  (t,  <J>),  and  then,  Jx  and  Jz  components  for  the  FD-TD  analysis  of  step  2. 
However,  for  the  sleeve-fitting  aperture,  only  the  coefficient  was  used 
to  form  J  (t,  <j>)  and  then,  for  the  FD-TD  analysis.  This  is  because 
|J^|  _<  0.05  | J t |  at  the  sleeve  fitting,  and  could  thus  be  neglected  without 
significant  error. 
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.  A-6.  GEOMETRY  AND  COORDINATE  SYSTEMS  FOR  MOM/FD-TD 
HYBRID  ANALYSIS  OF  MISSILE  GUIDANCE  SECTION, 
AXIAL  INCIDENCE 
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Fig.  A- 7.  MODE  COEFFICIENTS  FOR  APERTURE  EQUIVALENT  ELECTRIC  CURRENT,  HYBRID  MOM/FD-TD 
FORMULATION  OF  AXIAL  INCIDENCE  ILLUMINATION  OF  MISSILE  GUIDANCE  SECTION 


Integral  Equations  for  the  Finite  Wire  /  Narrow  Slot  Geometry 


This  appendix  summarizes  the  integral  equations  which  have  been  utilized 
to  obtain  the  numerical  results  shown  in  Section  7.0  of  this  report  for  the 
wire  /  slot  geometry  of  Figure  36.  The  complete  set  of  integro-differential 
equations  for  the  general  problem  of  a  finite-length  wire  excited  through  a 
narrow  slot  in  a  conducting  screen  appears  in  detail  in  the  literature  in 
References  11,  27,  and  28.  These  references  also  provide  numerical  results 
for  coupling  between  the  slot  and  the  wire  based  on  the  location  and  orienta¬ 
tion  of  the  wire  with  respect  to  the  slot. 

For  the  narrow  slot  /  finite  wire  geometry,  we  have  the  coupled  integral 
equations  [27]: 

'i'(x,)C{y')  K$s(x,  x'  ,  y' )  dx'dy' 

+  Jx!^  (X>  =  -3  f  Hx  (x)  (B-la) 
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e‘jkRww  g-jkR1 

(s>  s<)  = 


ww 


Rww  ’  [a2  + 


R1  .  [4z|  ♦  (S-S')2]1* 


3A 


3z  I  '  '2  '  cos6  '  2c  '  ft  / 
on  slot  s'=- 

(*,  y»0,  2-0)  2 


I  (s’)  G$w  (x,  s')  ds ’ 


ik  1  p”jkd 
GSW  (x’  s‘)  =  "  [X+  ~— 


d  =  {[x-(x„  +  s'  cos  a)]2  +  [y_  +  s’  cosB]2  +  [z  l2}*5 


3Fx  £  +  7 

'  =  zc  ST  /  2  m(x,)  {s>  x,)  dx’ 


3z 


on  wire 


JWS 


x'=-^ 


(B-2a) 

(B-2b) 

(B-3a) 

(B-3b) 

(B-3c) 

(B-4a) 

(B-4b) 

(B-4c) 

(B-5a) 


3WS 


( S ,  X '  )  =  -  [4£-  +  -y] 
u 


1 


-jkD 


(B-5b) 
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(B-5c) 


0  =  {[(xc  +  s  cos  a)  -  x']2  +  [y£  +  s  cos$]2  +  [z^2}1* 

The  integral  equations  (B-la  and  b)  are  solved  by  the  method-of-ewwhts 
numerical  approach  [13],  [27]. 
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